#####################################
# Install/load packages and functions
#####################################

# remove objects
rm(list=ls())
# detach all libraries
detachAllPackages <- function() {
  basic.packages <- c("package:stats", "package:graphics", "package:grDevices", "package:utils", "package:datasets", "package:methods", "package:base")
  package.list <- search()[ifelse(unlist(gregexpr("package:", search()))==1, TRUE, FALSE)]
  package.list <- setdiff(package.list, basic.packages)
  if (length(package.list)>0)  for (package in package.list) detach(package,  character.only=TRUE)
}
detachAllPackages()

# load libraries
pkgTest <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[,  "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg,  dependencies = TRUE)
  sapply(pkg,  require,  character.only = TRUE)
}

# here is where you load any necessary packages
# ex: stringr
# lapply(c("stringr"),  pkgTest)

lapply(c("car", "cregg", "data.table", "estimatr", "lmtest", "memisc", 
         "multcomp", "multiwayvcov", "pBrackets", "stringr", "texreg"),  pkgTest)

# set working directory
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# set seed in honor of T-Swift
set.seed(1989)

###################
# load and fix data 
###################

# load data without first two rows of unneccessary variable names
data <- fread("mummolo_data.csv", header = TRUE, stringsAsFactors = FALSE)
data <- data[-c(1:2),]

# The coding scheme used by Mummolo requires the researcher to identify which of
# six affected publics the respondent might belong to--women, smokers, senior
# citizens, college students, uninsured persons, healthcare workers, or 
# those trying to lose weight.  The following code chunks create binary indicators
# for whether each respondent belongs to each affected public.

data$ap_women <- ifelse(data$Q4=="Female", 1, NA)
data$ap_women <- ifelse(data$Q4=="Male" | data$Q4=="Other (Please specify)", 0, data$ap_women)

data$ap_smokers <- ifelse(data$Q20=="I never smoke" | data$Q20=="One cigarette", 0, NA)
data$ap_smokers <- ifelse(data$Q20=="One pack" | data$Q20=="Two to four packs" | data$Q20=="Five to ten packs" | data$Q20=="More than ten packs", 1, data$ap_smokers)

data$ap_weight <- ifelse(data$Q19=="Yes", 1, NA)
data$ap_weight <- ifelse(data$Q19=="No", 0, data$ap_weight)

data$ap_seniors <- ifelse(data$Q2<1966, 1, NA)
data$ap_seniors <- ifelse(data$Q2>=1966, 0, data$ap_seniors)

data$ap_students <- ifelse(data$Q18=="Yes, I'm a full-time student" | data$Q18=="Yes, I'm a part-time student", 1, NA)
data$ap_students <- ifelse(data$Q18=="No, I am not currently a student", 0, data$ap_students)

data$ap_health <- ifelse(data$Q16=="Yes" & data$Q17=="No", 0, NA)
data$ap_health <- ifelse(data$Q16=="No" | data$Q17=="Yes", 1, data$ap_health)

# Following Mummolo, we code source friendliness as the correspondence of PID
# and the perceived partisan slant of the news organization featured.  Here,
# we recode respondents' PID to an ordinal scale

data$pid7 <- ifelse(data$Q13=="Strong Democrat", 7, NA)
data$pid7 <- ifelse(data$Q13=="Not a very strong Democrat", 6, data$pid7)
data$pid7 <- ifelse(data$Q15=="Democratic", 5, data$pid7)
data$pid7 <- ifelse(data$Q15=="Neither" | data$Q12=="Other", 4, data$pid7)
data$pid7 <- ifelse(data$Q15=="Republican", 3, data$pid7)
data$pid7 <- ifelse(data$Q14=="Not a very strong Republican", 2, data$pid7)
data$pid7 <- ifelse(data$Q14=="Strong Republican", 1, data$pid7)

# Here, we code which experimental arm each respondent was in and which type of
# task these consequently encountered in their 13th task
data$exp_arm <- ifelse(data$FL_8_DO=="FL_15", "Forced", "Abstention")
data$final_task <- ifelse(data$FL_20_DO!="" & data$FL_50_DO=="", "Forced", NA)
data$final_task <- ifelse(data$FL_20_DO=="" & data$FL_50_DO!="", "Abstention", data$final_task)

# For each of the following code chunks, we determine for each task which conjoint
# attributes correspond with source and headline and organize that information
# for each profile, as well as recover the respondent choice
data$headline1_A <- ifelse(data$`F-1-1`=="Headline", data$`F-1-1-1`, data$`F-1-1-2`)
data$source1_A <- ifelse(data$`F-1-1`=="Source", data$`F-1-1-1`, data$`F-1-1-2`)
data$headline1_B <- ifelse(data$`F-1-1`=="Headline", data$`F-1-2-1`, data$`F-1-2-2`)
data$source1_B <- ifelse(data$`F-1-1`=="Source", data$`F-1-2-1`, data$`F-1-2-2`)
data$choice1 <- ifelse(data$exp_arm=="Forced", data$Q23, data$Q52)

data$headline2_A <- ifelse(data$`F-2-1`=="Headline", data$`F-2-1-1`, data$`F-2-1-2`)
data$source2_A <- ifelse(data$`F-2-1`=="Source", data$`F-2-1-1`, data$`F-2-1-2`)
data$headline2_B <- ifelse(data$`F-2-1`=="Headline", data$`F-2-2-1`, data$`F-2-2-2`)
data$source2_B <- ifelse(data$`F-2-1`=="Source", data$`F-2-2-1`, data$`F-2-2-2`)
data$choice2 <- ifelse(data$exp_arm=="Forced", data$Q25, data$Q54)

data$headline3_A <- ifelse(data$`F-3-1`=="Headline", data$`F-3-1-1`, data$`F-3-1-2`)
data$source3_A <- ifelse(data$`F-3-1`=="Source", data$`F-3-1-1`, data$`F-3-1-2`)
data$headline3_B <- ifelse(data$`F-3-1`=="Headline", data$`F-3-2-1`, data$`F-3-2-2`)
data$source3_B <- ifelse(data$`F-3-1`=="Source", data$`F-3-2-1`, data$`F-3-2-2`)
data$choice3 <- ifelse(data$exp_arm=="Forced", data$Q29, data$Q56)

data$headline4_A <- ifelse(data$`F-4-1`=="Headline", data$`F-4-1-1`, data$`F-4-1-2`)
data$source4_A <- ifelse(data$`F-4-1`=="Source", data$`F-4-1-1`, data$`F-4-1-2`)
data$headline4_B <- ifelse(data$`F-4-1`=="Headline", data$`F-4-2-1`, data$`F-4-2-2`)
data$source4_B <- ifelse(data$`F-4-1`=="Source", data$`F-4-2-1`, data$`F-4-2-2`)
data$choice4 <- ifelse(data$exp_arm=="Forced", data$Q31, data$Q58)

data$headline5_A <- ifelse(data$`F-5-1`=="Headline", data$`F-5-1-1`, data$`F-5-1-2`)
data$source5_A <- ifelse(data$`F-5-1`=="Source", data$`F-5-1-1`, data$`F-5-1-2`)
data$headline5_B <- ifelse(data$`F-5-1`=="Headline", data$`F-5-2-1`, data$`F-5-2-2`)
data$source5_B <- ifelse(data$`F-5-1`=="Source", data$`F-5-2-1`, data$`F-5-2-2`)
data$choice5 <- ifelse(data$exp_arm=="Forced", data$Q33, data$Q60)

data$headline6_A <- ifelse(data$`F-6-1`=="Headline", data$`F-6-1-1`, data$`F-6-1-2`)
data$source6_A <- ifelse(data$`F-6-1`=="Source", data$`F-6-1-1`, data$`F-6-1-2`)
data$headline6_B <- ifelse(data$`F-6-1`=="Headline", data$`F-6-2-1`, data$`F-6-2-2`)
data$source6_B <- ifelse(data$`F-6-1`=="Source", data$`F-6-2-1`, data$`F-6-2-2`)
data$choice6 <- ifelse(data$exp_arm=="Forced", data$Q35, data$Q62)

data$headline7_A <- ifelse(data$`F-7-1`=="Headline", data$`F-7-1-1`, data$`F-7-1-2`)
data$source7_A <- ifelse(data$`F-7-1`=="Source", data$`F-7-1-1`, data$`F-7-1-2`)
data$headline7_B <- ifelse(data$`F-7-1`=="Headline", data$`F-7-2-1`, data$`F-7-2-2`)
data$source7_B <- ifelse(data$`F-7-1`=="Source", data$`F-7-2-1`, data$`F-7-2-2`)
data$choice7 <- ifelse(data$exp_arm=="Forced", data$Q37, data$Q64)

data$headline8_A <- ifelse(data$`F-8-1`=="Headline", data$`F-8-1-1`, data$`F-8-1-2`)
data$source8_A <- ifelse(data$`F-8-1`=="Source", data$`F-8-1-1`, data$`F-8-1-2`)
data$headline8_B <- ifelse(data$`F-8-1`=="Headline", data$`F-8-2-1`, data$`F-8-2-2`)
data$source8_B <- ifelse(data$`F-8-1`=="Source", data$`F-8-2-1`, data$`F-8-2-2`)
data$choice8 <- ifelse(data$exp_arm=="Forced", data$Q39, data$Q66)

data$headline9_A <- ifelse(data$`F-9-1`=="Headline", data$`F-9-1-1`, data$`F-9-1-2`)
data$source9_A <- ifelse(data$`F-9-1`=="Source", data$`F-9-1-1`, data$`F-9-1-2`)
data$headline9_B <- ifelse(data$`F-9-1`=="Headline", data$`F-9-2-1`, data$`F-9-2-2`)
data$source9_B <- ifelse(data$`F-9-1`=="Source", data$`F-9-2-1`, data$`F-9-2-2`)
data$choice9 <- ifelse(data$exp_arm=="Forced", data$Q41, data$Q68)

data$headline10_A <- ifelse(data$`F-10-1`=="Headline", data$`F-10-1-1`, data$`F-10-1-2`)
data$source10_A <- ifelse(data$`F-10-1`=="Source", data$`F-10-1-1`, data$`F-10-1-2`)
data$headline10_B <- ifelse(data$`F-10-1`=="Headline", data$`F-10-2-1`, data$`F-10-2-2`)
data$source10_B <- ifelse(data$`F-10-1`=="Source", data$`F-10-2-1`, data$`F-10-2-2`)
data$choice10 <- ifelse(data$exp_arm=="Forced", data$Q43, data$Q70)

data$headline11_A <- ifelse(data$`F-11-1`=="Headline", data$`F-11-1-1`, data$`F-11-1-2`)
data$source11_A <- ifelse(data$`F-11-1`=="Source", data$`F-11-1-1`, data$`F-11-1-2`)
data$headline11_B <- ifelse(data$`F-11-1`=="Headline", data$`F-11-2-1`, data$`F-11-2-2`)
data$source11_B <- ifelse(data$`F-11-1`=="Source", data$`F-11-2-1`, data$`F-11-2-2`)
data$choice11 <- ifelse(data$exp_arm=="Forced", data$Q45, data$Q72)

data$headline12_A <- ifelse(data$`F-12-1`=="Headline", data$`F-12-1-1`, data$`F-12-1-2`)
data$source12_A <- ifelse(data$`F-12-1`=="Source", data$`F-12-1-1`, data$`F-12-1-2`)
data$headline12_B <- ifelse(data$`F-12-1`=="Headline", data$`F-12-2-1`, data$`F-12-2-2`)
data$source12_B <- ifelse(data$`F-12-1`=="Source", data$`F-12-2-1`, data$`F-12-2-2`)
data$choice12 <- ifelse(data$exp_arm=="Forced", data$Q27, data$Q74)

# Before reorganizing the conjoint attribute-levels and choices, we save the 
# respondent-level information
demo_data <- subset(data, select = c(ResponseId, ap_women, ap_smokers, ap_health, ap_weight, ap_seniors, ap_students,
                                     pid7, exp_arm))

# The following code chunks make a new long dataframe for each profile, which we
# then rbind together as a long dataframe.
long_data1_A <- subset(data, select = c(ResponseId, headline1_A, source1_A, choice1))
long_data1_A$task <- 1
long_data1_A$profile <- 1
colnames(long_data1_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data1_B <- subset(data, select = c(ResponseId, headline1_B, source1_B, choice1))
long_data1_B$task <- 1
long_data1_B$profile <- 2
colnames(long_data1_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data2_A <- subset(data, select = c(ResponseId, headline2_A, source2_A, choice2))
long_data2_A$task <- 2
long_data2_A$profile <- 1
colnames(long_data2_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data2_B <- subset(data, select = c(ResponseId, headline2_B, source2_B, choice2))
long_data2_B$task <- 2
long_data2_B$profile <- 2
colnames(long_data2_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data3_A <- subset(data, select = c(ResponseId, headline3_A, source3_A, choice3))
long_data3_A$task <- 3
long_data3_A$profile <- 1
colnames(long_data3_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data3_B <- subset(data, select = c(ResponseId, headline3_B, source3_B, choice3))
long_data3_B$task <- 3
long_data3_B$profile <- 2
colnames(long_data3_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data4_A <- subset(data, select = c(ResponseId, headline4_A, source4_A, choice4))
long_data4_A$task <- 4
long_data4_A$profile <- 1
colnames(long_data4_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data4_B <- subset(data, select = c(ResponseId, headline4_B, source4_B, choice4))
long_data4_B$task <- 4
long_data4_B$profile <- 2
colnames(long_data4_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data5_A <- subset(data, select = c(ResponseId, headline5_A, source5_A, choice5))
long_data5_A$task <- 5
long_data5_A$profile <- 1
colnames(long_data5_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data5_B <- subset(data, select = c(ResponseId, headline5_B, source5_B, choice5))
long_data5_B$task <- 5
long_data5_B$profile <- 2
colnames(long_data5_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data6_A <- subset(data, select = c(ResponseId, headline6_A, source6_A, choice6))
long_data6_A$task <- 6
long_data6_A$profile <- 1
colnames(long_data6_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data6_B <- subset(data, select = c(ResponseId, headline6_B, source6_B, choice6))
long_data6_B$task <- 6
long_data6_B$profile <- 2
colnames(long_data6_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data7_A <- subset(data, select = c(ResponseId, headline7_A, source7_A, choice7))
long_data7_A$task <- 7
long_data7_A$profile <- 1
colnames(long_data7_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data7_B <- subset(data, select = c(ResponseId, headline7_B, source7_B, choice7))
long_data7_B$task <- 7
long_data7_B$profile <- 2
colnames(long_data7_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data8_A <- subset(data, select = c(ResponseId, headline8_A, source8_A, choice8))
long_data8_A$task <- 8
long_data8_A$profile <- 1
colnames(long_data8_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data8_B <- subset(data, select = c(ResponseId, headline8_B, source8_B, choice8))
long_data8_B$task <- 8
long_data8_B$profile <- 2
colnames(long_data8_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data9_A <- subset(data, select = c(ResponseId, headline9_A, source9_A, choice9))
long_data9_A$task <- 9
long_data9_A$profile <- 1
colnames(long_data9_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data9_B <- subset(data, select = c(ResponseId, headline9_B, source9_B, choice9))
long_data9_B$task <- 9
long_data9_B$profile <- 2
colnames(long_data9_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data10_A <- subset(data, select = c(ResponseId, headline10_A, source10_A, choice10))
long_data10_A$task <- 10
long_data10_A$profile <- 1
colnames(long_data10_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data10_B <- subset(data, select = c(ResponseId, headline10_B, source10_B, choice10))
long_data10_B$task <- 10
long_data10_B$profile <- 2
colnames(long_data10_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data11_A <- subset(data, select = c(ResponseId, headline11_A, source11_A, choice11))
long_data11_A$task <- 11
long_data11_A$profile <- 1
colnames(long_data11_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data11_B <- subset(data, select = c(ResponseId, headline11_B, source11_B, choice11))
long_data11_B$task <- 11
long_data11_B$profile <- 2
colnames(long_data11_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data12_A <- subset(data, select = c(ResponseId, headline12_A, source12_A, choice12))
long_data12_A$task <- 12
long_data12_A$profile <- 1
colnames(long_data12_A) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data12_B <- subset(data, select = c(ResponseId, headline12_B, source12_B, choice12))
long_data12_B$task <- 12
long_data12_B$profile <- 2
colnames(long_data12_B) <- c("ResponseId", "headline", "source", "choice", "task", "profile")

long_data <- rbind(long_data1_A, long_data1_B, long_data2_A, long_data2_B, long_data3_A, long_data3_B, long_data4_A, long_data4_B,
                   long_data5_A, long_data5_B, long_data6_A, long_data6_B, long_data7_A, long_data7_B, long_data8_A, long_data8_B,
                   long_data9_A, long_data9_B, long_data10_A, long_data10_B, long_data11_A, long_data11_B, long_data12_A, long_data12_B)

# Merging the respondent-level data back into the new long dataframe.
full_data <- merge(demo_data, long_data, by = "ResponseId")

# In the following code chunks, we use information about the profile's source and
# headline, as well as the respondents' PID and membership in the affected
# publics, to determine whether the source is friendly, unfriendly, or neutral and
# whether the headline is relevant or irrelevant
full_data$source_friendliness <- ifelse(full_data$source=="MSNBC" & full_data$pid7>4, "Friendly", NA)
full_data$source_friendliness <- ifelse(full_data$source=="MSNBC" & full_data$pid7<4, "Unfriendly", full_data$source_friendliness)
full_data$source_friendliness <- ifelse(full_data$source=="MSNBC" & full_data$pid7==4, "Neutral", full_data$source_friendliness)
full_data$source_friendliness <- ifelse(full_data$source=="Fox News" & full_data$pid7<4, "Friendly", full_data$source_friendliness)
full_data$source_friendliness <- ifelse(full_data$source=="Fox News" & full_data$pid7>4, "Unfriendly", full_data$source_friendliness)
full_data$source_friendliness <- ifelse(full_data$source=="Fox News" & full_data$pid7==4, "Neutral", full_data$source_friendliness)
full_data$source_friendliness <- ifelse(full_data$source=="USA Today" & full_data$pid7<4, "Neutral", full_data$source_friendliness)
full_data$source_friendliness <- ifelse(full_data$source=="USA Today" & full_data$pid7>4, "Neutral", full_data$source_friendliness)
full_data$source_friendliness <- ifelse(full_data$source=="USA Today" & full_data$pid7==4, "Neutral", full_data$source_friendliness)

full_data$topic_relevance <- ifelse(full_data$headline=="2021 MLB Baseball Season Guide", "Irrelevant", NA)
full_data$topic_relevance <- ifelse(full_data$headline=="2021 MLB Baseball Season Outlook", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="2021 MLB Baseball Season Preview", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Abortion restrictions reinstated by federal court" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Abortion restrictions reinstated by federal court" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Bill to guarantee equal pay for women blocked in Senate" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Bill to guarantee equal pay for women blocked in Senate" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Bill to mandate equal pay for women fails in Senate vote" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Bill to mandate equal pay for women fails in Senate vote" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Celebrities' worst dating stories", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Celebrity dating fails", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Celebrity dating nightmares", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Federal court puts abortion restrictions back in place" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Federal court puts abortion restrictions back in place" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Federal court reinstates abortion restrictions" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Federal court reinstates abortion restrictions" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Congress considering cuts to Social Security" & full_data$ap_seniors==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Congress considering cuts to Social Security" & full_data$ap_seniors==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Congress debates cuts to Social Security" & full_data$ap_seniors==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Congress debates cuts to Social Security" & full_data$ap_seniors==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Congress weighs cuts to Social Security" & full_data$ap_seniors==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Congress weighs cuts to Social Security" & full_data$ap_seniors==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="How to drop all that student loan debt" & full_data$ap_students==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="How to drop all that student loan debt" & full_data$ap_students==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Here are the new members of Joe Biden's Cabinet", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Introducing Joe Biden's new Cabinet members", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Meet the new members of Joe Biden's Cabinet", "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Obamacare enrollment figures met with doubt" & full_data$ap_health==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Obamacare enrollment figures met with doubt" & full_data$ap_health==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Obamacare enrollment numbers called into question" & full_data$ap_health==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Obamacare enrollment numbers called into question" & full_data$ap_health==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Some question Obamacare enrollment figures" & full_data$ap_health==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Some question Obamacare enrollment figures" & full_data$ap_health==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Senate votes against bill that would ensure equal pay for women" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Senate votes against bill that would ensure equal pay for women" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Smokers who quit may cut heart risk faster than had been thought, study finds" & full_data$ap_smokers==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Smokers who quit may cut heart risk faster than had been thought, study finds" & full_data$ap_smokers==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study finds potential cause of breast cancer" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study finds potential cause of breast cancer" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study finds those who quit smoking repair heart damage fast than expected" & full_data$ap_smokers==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study finds those who quit smoking repair heart damage fast than expected" & full_data$ap_smokers==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study identifies possible breast cancer cause" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study identifies possible breast cancer cause" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study points to new cause of breast cancer" & full_data$ap_women==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study points to new cause of breast cancer" & full_data$ap_women==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study: quitting smoking repairs heart damage more quickly than was thought" & full_data$ap_smokers==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Study: quitting smoking repairs heart damage more quickly than was thought" & full_data$ap_smokers==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="The best ways to get rid of student loan debt" & full_data$ap_students==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="The best ways to get rid of student loan debt" & full_data$ap_students==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Tips for getting rid of student loan debt" & full_data$ap_students==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Tips for getting rid of student loan debt" & full_data$ap_students==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Weight-loss tips and success stories" & full_data$ap_weight==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Weight-loss tips and success stories" & full_data$ap_weight==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Weight-loss tips that make a difference" & full_data$ap_weight==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Weight-loss tips that make a difference" & full_data$ap_weight==0, "Irrelevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Weight-loss tips that work" & full_data$ap_weight==1, "Relevant", full_data$topic_relevance)
full_data$topic_relevance <- ifelse(full_data$headline=="Weight-loss tips that work" & full_data$ap_weight==0, "Irrelevant", full_data$topic_relevance)

# Here, we recode respondents' choices as a binary variable.
full_data$choice_bin <- ifelse(full_data$choice=="News Selection A" & full_data$profile==1, 1, NA)
full_data$choice_bin <- ifelse(full_data$choice=="News Selection A" & full_data$profile==2, 0, full_data$choice_bin)
full_data$choice_bin <- ifelse(full_data$choice=="News Selection B" & full_data$profile==2, 1, full_data$choice_bin)
full_data$choice_bin <- ifelse(full_data$choice=="News Selection B" & full_data$profile==1, 0, full_data$choice_bin)
full_data$choice_bin <- ifelse(full_data$choice=="Neither", 0, full_data$choice_bin)

# Here, we specify the levels and order of those levels for our key factor
# variables
full_data$source_friendliness <- factor(full_data$source_friendliness, levels = c("Neutral", "Friendly", "Unfriendly"))
full_data$topic_relevance <- factor(full_data$topic_relevance, levels = c("Irrelevant", "Relevant"))
full_data$exp_arm <- factor(full_data$exp_arm, levels = c("Forced", "Abstention"))

# Mummolo retained only partisans for analysis, so all who identified as true 
# independents need to be dropped

full_data_part <- full_data[which(full_data$pid7!=4),]

#############
# FIGURE SM.5
#############

# In each of the following code chunks, we estimate a model for one of the 
# six issue areas meant to cue an attentive public, then obtain the
# cluster-robust variance-covariance matrix (clustered on respondent) and
# reassign that matrix to the model object and save the standard errors and 
# p-values

women_int_reg <- lm(choice_bin ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_women==1),])
summary(women_int_reg)
women_int_reg_vcovCL <- cluster.vcov(women_int_reg, full_data_part$ResponseId[which(full_data_part$ap_women==1)])
coeftest(women_int_reg, women_int_reg_vcovCL)
women_int_reg <-withVCov(women_int_reg, vcov = women_int_reg_vcovCL)
women_int_reg_ses <- coeftest(women_int_reg, women_int_reg_vcovCL)[,2]
women_int_reg_pvals <- coeftest(women_int_reg, women_int_reg_vcovCL)[,4]

smokers_int_reg <- lm(choice_bin ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_smokers==1),])
summary(smokers_int_reg)
smokers_int_reg_vcovCL <- cluster.vcov(smokers_int_reg, full_data_part$ResponseId[which(full_data_part$ap_smokers==1)])
coeftest(smokers_int_reg, smokers_int_reg_vcovCL)
smokers_int_reg <-withVCov(smokers_int_reg, vcov = smokers_int_reg_vcovCL)
smokers_int_reg_ses <- coeftest(smokers_int_reg, smokers_int_reg_vcovCL)[,2]
smokers_int_reg_pvals <- coeftest(smokers_int_reg, smokers_int_reg_vcovCL)[,4]

seniors_int_reg <- lm(choice_bin ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_seniors==1),])
summary(seniors_int_reg)
seniors_int_reg_vcovCL <- cluster.vcov(seniors_int_reg, full_data_part$ResponseId[which(full_data_part$ap_seniors==1)])
coeftest(seniors_int_reg, seniors_int_reg_vcovCL)
seniors_int_reg <-withVCov(seniors_int_reg, vcov = seniors_int_reg_vcovCL)
seniors_int_reg_ses <- coeftest(seniors_int_reg, seniors_int_reg_vcovCL)[,2]
seniors_int_reg_pvals <- coeftest(seniors_int_reg, seniors_int_reg_vcovCL)[,4]

students_int_reg <- lm(choice_bin ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_students==1),])
summary(students_int_reg)
students_int_reg_vcovCL <- cluster.vcov(students_int_reg, full_data_part$ResponseId[which(full_data_part$ap_students==1)])
coeftest(students_int_reg, students_int_reg_vcovCL)
students_int_reg <-withVCov(students_int_reg, vcov = students_int_reg_vcovCL)
students_int_reg_ses <- coeftest(students_int_reg, students_int_reg_vcovCL)[,2]
students_int_reg_pvals <- coeftest(students_int_reg, students_int_reg_vcovCL)[,4]

health_int_reg <- lm(choice_bin ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_health==1),])
summary(health_int_reg)
health_int_reg_vcovCL <- cluster.vcov(health_int_reg, full_data_part$ResponseId[which(full_data_part$ap_health==1)])
coeftest(health_int_reg, health_int_reg_vcovCL)
health_int_reg <-withVCov(health_int_reg, vcov = health_int_reg_vcovCL)
health_int_reg_ses <- coeftest(health_int_reg, health_int_reg_vcovCL)[,2]
health_int_reg_pvals <- coeftest(health_int_reg, health_int_reg_vcovCL)[,4]

weight_int_reg <- lm(choice_bin ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_weight==1),])
summary(weight_int_reg)
weight_int_reg_vcovCL <- cluster.vcov(weight_int_reg, full_data_part$ResponseId[which(full_data_part$ap_weight==1)])
coeftest(weight_int_reg, weight_int_reg_vcovCL)
weight_int_reg <-withVCov(weight_int_reg, vcov = weight_int_reg_vcovCL)
weight_int_reg_ses <- coeftest(weight_int_reg, weight_int_reg_vcovCL)[,2]
weight_int_reg_pvals <- coeftest(weight_int_reg, weight_int_reg_vcovCL)[,4]

# In the following pairs of code chunks, we estimate the AMCEs associated with
# those in the abstention and forced choice arms for each unique pairing of
# friendliness-topic relevance, then estimate the under/overestimation of the
# forced choice AMCEs relative to the abstention option AMCEs

women_amces <- confint(glht(women_int_reg, 
                            linfct = c("source_friendlinessUnfriendly = 0",
                                       "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                       "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                       "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                       "topic_relevanceRelevant = 0",
                                       "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                       "source_friendlinessFriendly = 0",
                                       "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                       "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                       "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                       calpha = univariate_calpha())$confint

women_amce_diffs <- confint(glht(women_int_reg, 
                                 linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                            "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                            "topic_relevanceRelevant:exp_armAbstention = 0",
                                            "source_friendlinessFriendly:exp_armAbstention = 0",
                                            "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                            calpha = 2.94)$confint

smokers_amces <- confint(glht(smokers_int_reg, 
                              linfct = c("source_friendlinessUnfriendly = 0",
                                         "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                         "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                         "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                         "topic_relevanceRelevant = 0",
                                         "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                         "source_friendlinessFriendly = 0",
                                         "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                         "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                         "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                         calpha = univariate_calpha())$confint

smokers_amce_diffs <- confint(glht(smokers_int_reg, 
                                   linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                              "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                              "topic_relevanceRelevant:exp_armAbstention = 0",
                                              "source_friendlinessFriendly:exp_armAbstention = 0",
                                              "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                              calpha = 2.94)$confint

seniors_amces <- confint(glht(seniors_int_reg, 
                              linfct = c("source_friendlinessUnfriendly = 0",
                                         "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                         "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                         "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                         "topic_relevanceRelevant = 0",
                                         "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                         "source_friendlinessFriendly = 0",
                                         "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                         "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                         "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                         calpha = univariate_calpha())$confint

seniors_amce_diffs <- confint(glht(seniors_int_reg, 
                                   linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                              "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                              "topic_relevanceRelevant:exp_armAbstention = 0",
                                              "source_friendlinessFriendly:exp_armAbstention = 0",
                                              "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                              calpha = 2.94)$confint

students_amces <- confint(glht(students_int_reg, 
                               linfct = c("source_friendlinessUnfriendly = 0",
                                          "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                          "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                          "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                          "topic_relevanceRelevant = 0",
                                          "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                          "source_friendlinessFriendly = 0",
                                          "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                          "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                          "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                          calpha = univariate_calpha())$confint

students_amce_diffs <- confint(glht(students_int_reg, 
                                    linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                               "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                               "topic_relevanceRelevant:exp_armAbstention = 0",
                                               "source_friendlinessFriendly:exp_armAbstention = 0",
                                               "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                               calpha = 2.94)$confint

health_amces <- confint(glht(health_int_reg, 
                             linfct = c("source_friendlinessUnfriendly = 0",
                                        "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                        "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                        "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                        "topic_relevanceRelevant = 0",
                                        "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                        "source_friendlinessFriendly = 0",
                                        "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                        "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                        "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                        calpha = univariate_calpha())$confint

health_amce_diffs <- confint(glht(health_int_reg, 
                                  linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                             "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                             "topic_relevanceRelevant:exp_armAbstention = 0",
                                             "source_friendlinessFriendly:exp_armAbstention = 0",
                                             "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                             calpha = 2.94)$confint

weight_amces <- confint(glht(weight_int_reg, 
                             linfct = c("source_friendlinessUnfriendly = 0",
                                        "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                        "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                        "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                        "topic_relevanceRelevant = 0",
                                        "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                        "source_friendlinessFriendly = 0",
                                        "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                        "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                        "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                        calpha = univariate_calpha())$confint

weight_amce_diffs <- confint(glht(weight_int_reg, 
                                  linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                             "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                             "topic_relevanceRelevant:exp_armAbstention = 0",
                                             "source_friendlinessFriendly:exp_armAbstention = 0",
                                             "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                             calpha = 2.94)$confint

# NOTE: In the following code that creates the figure, some of the differences
# in AMCEs from the forced choice to abstention arms are multiplied by
# negative 1; when this is the case, it is to keep consistent the meaning of
# the values presented in the right pane of the figure difference in the magnitudes
# of the two AMCEs irrespective of signs)

pdf(file = "../figures/figure_SM5.pdf", family = "Times", height = 18, width=13)
layout(matrix(c(1,2), nrow=2, byrow=TRUE), heights = c(0.85,0.15))
par(mar=c(6.1,12,5.75,.5))
plot(x=c(weight_amces[10,1],
         weight_amces[9,1],
         weight_amces[8,1],
         weight_amces[7,1],
         weight_amces[6,1],
         weight_amces[5,1],
         0,
         0,
         weight_amces[4,1],
         weight_amces[3,1],
         weight_amces[2,1],
         weight_amces[1,1],
         health_amces[10,1],
         health_amces[9,1],
         health_amces[8,1],
         health_amces[7,1],
         health_amces[6,1],
         health_amces[5,1],
         0,
         0,
         health_amces[4,1],
         health_amces[3,1],
         health_amces[2,1],
         health_amces[1,1],
         students_amces[10,1],
         students_amces[9,1],
         students_amces[8,1],
         students_amces[7,1],
         students_amces[6,1],
         students_amces[5,1],
         0,
         0,
         students_amces[4,1],
         students_amces[3,1],
         students_amces[2,1],
         students_amces[1,1],
         seniors_amces[10,1],
         seniors_amces[9,1],
         seniors_amces[8,1],
         seniors_amces[7,1],
         seniors_amces[6,1],
         seniors_amces[5,1],
         0,
         0,
         seniors_amces[4,1],
         seniors_amces[3,1],
         seniors_amces[2,1],
         seniors_amces[1,1],
         smokers_amces[10,1],
         smokers_amces[9,1],
         smokers_amces[8,1],
         smokers_amces[7,1],
         smokers_amces[6,1],
         smokers_amces[5,1],
         0,
         0,
         smokers_amces[4,1],
         smokers_amces[3,1],
         smokers_amces[2,1],
         smokers_amces[1,1],
         women_amces[10,1],
         women_amces[9,1],
         women_amces[8,1],
         women_amces[7,1],
         women_amces[6,1],
         women_amces[5,1],
         0,
         0,
         women_amces[4,1],
         women_amces[3,1],
         women_amces[2,1],
         women_amces[1,1],
         -weight_amce_diffs[5,1]+0.85,
         -weight_amce_diffs[4,1]+0.85,
         -weight_amce_diffs[3,1]+0.85,
         0+0.85,
         -weight_amce_diffs[2,1]+0.85,
         -weight_amce_diffs[1,1]+0.85,
         -health_amce_diffs[5,1]+0.85,
         -health_amce_diffs[4,1]+0.85,
         -health_amce_diffs[3,1]+0.85,
         0+0.85,
         -health_amce_diffs[2,1]+0.85,
         health_amce_diffs[1,1]+0.85,
         -students_amce_diffs[5,1]+0.85,
         -students_amce_diffs[4,1]+0.85,
         -students_amce_diffs[3,1]+0.85,
         0+0.85,
         -students_amce_diffs[2,1]+0.85,
         students_amce_diffs[1,1]+0.85,
         -seniors_amce_diffs[5,1]+0.85,
         seniors_amce_diffs[4,1]+0.85,
         -seniors_amce_diffs[3,1]+0.85,
         0+0.85,
         -seniors_amce_diffs[2,1]+0.85,
         seniors_amce_diffs[1,1]+0.85,
         -smokers_amce_diffs[5,1]+0.85,
         -smokers_amce_diffs[4,1]+0.85,
         -smokers_amce_diffs[3,1]+0.85,
         0+0.85,
         -smokers_amce_diffs[2,1]+0.85,
         smokers_amce_diffs[1,1]+0.85,
         -women_amce_diffs[5,1]+0.85,
         -women_amce_diffs[4,1]+0.85,
         -women_amce_diffs[3,1]+0.85,
         0+0.85,
         -women_amce_diffs[2,1]+0.85,
         women_amce_diffs[1,1]+0.85), 
     y=c(1.3, 1.7,
         2.3, 2.7,
         3.3, 3.7,
         4.3, 4.7,
         5.3, 5.7,
         6.3, 6.7,
         8.3, 8.7,
         9.3, 9.7,
         10.3, 10.7,
         11.3, 11.7,
         12.3, 12.7,
         13.3, 13.7,
         15.3, 15.7,
         16.3, 16.7,
         17.3, 17.7,
         18.3, 18.7,
         19.3, 19.7,
         20.3, 20.7,
         22.3, 22.7,
         23.3, 23.7,
         24.3, 24.7,
         25.3, 25.7,
         26.3, 26.7,
         27.3, 27.7,
         29.3, 29.7,
         30.3, 30.7,
         31.3, 31.7,
         32.3, 32.7,
         33.3, 33.7,
         34.3, 34.7,
         36.3, 36.7,
         37.3, 37.7,
         38.3, 38.7,
         39.3, 39.7,
         40.3, 40.7,
         41.3, 41.7,
         seq(1.5,6.5,1),
         seq(8.5,13.5,1),
         seq(15.5,20.5,1),
         seq(22.5,27.5,1),
         seq(29.5,34.5,1),
         seq(36.5,41.5,1)),
     xlim=c(-0.25,1.15),
     xaxt="n",
     ylim=c(1, 42),
     pch=c(rep(c(17,19),36), rep(15, 36)),
     col=c(rep(c("gray50","black"),36), rep("black", 36)),
     tck=-.02,
     cex.axis=0.9,
     cex=1.25,
     ylab="",
     yaxt="n",
     xlab="",
     axes = FALSE,
     panel.first = c(abline(v=0,lwd=2, col="gray70",lty=2),
                     abline(v=0.85,lwd=2, col="gray70",lty=2)))
segments(col=c(rep(c("gray50", "black"),36), rep("black", 36)),
         x0=c(weight_amces[10,2],
              weight_amces[9,2],
              weight_amces[8,2],
              weight_amces[7,2],
              weight_amces[6,2],
              weight_amces[5,2],
              0,
              0,
              weight_amces[4,2],
              weight_amces[3,2],
              weight_amces[2,2],
              weight_amces[1,2],
              health_amces[10,2],
              health_amces[9,2],
              health_amces[8,2],
              health_amces[7,2],
              health_amces[6,2],
              health_amces[5,2],
              0,
              0,
              health_amces[4,2],
              health_amces[3,2],
              health_amces[2,2],
              health_amces[1,2],
              students_amces[10,2],
              students_amces[9,2],
              students_amces[8,2],
              students_amces[7,2],
              students_amces[6,2],
              students_amces[5,2],
              0,
              0,
              students_amces[4,2],
              students_amces[3,2],
              students_amces[2,2],
              students_amces[1,2],
              seniors_amces[10,2],
              seniors_amces[9,2],
              seniors_amces[8,2],
              seniors_amces[7,2],
              seniors_amces[6,2],
              seniors_amces[5,2],
              0,
              0,
              seniors_amces[4,2],
              seniors_amces[3,2],
              seniors_amces[2,2],
              seniors_amces[1,2],
              smokers_amces[10,2],
              smokers_amces[9,2],
              smokers_amces[8,2],
              smokers_amces[7,2],
              smokers_amces[6,2],
              smokers_amces[5,2],
              0,
              0,
              smokers_amces[4,2],
              smokers_amces[3,2],
              smokers_amces[2,2],
              smokers_amces[1,2],
              women_amces[10,2],
              women_amces[9,2],
              women_amces[8,2],
              women_amces[7,2],
              women_amces[6,2],
              women_amces[5,2],
              0,
              0,
              women_amces[4,2],
              women_amces[3,2],
              women_amces[2,2],
              women_amces[1,2],
              -weight_amce_diffs[5,2]+0.85,
              -weight_amce_diffs[4,2]+0.85,
              -weight_amce_diffs[3,2]+0.85,
              0+0.85,
              -weight_amce_diffs[2,2]+0.85,
              -weight_amce_diffs[1,2]+0.85,
              -health_amce_diffs[5,2]+0.85,
              -health_amce_diffs[4,2]+0.85,
              -health_amce_diffs[3,2]+0.85,
              0+0.85,
              -health_amce_diffs[2,2]+0.85,
              health_amce_diffs[1,2]+0.85,
              -students_amce_diffs[5,2]+0.85,
              -students_amce_diffs[4,2]+0.85,
              -students_amce_diffs[3,2]+0.85,
              0+0.85,
              -students_amce_diffs[2,2]+0.85,
              students_amce_diffs[1,2]+0.85,
              -seniors_amce_diffs[5,2]+0.85,
              seniors_amce_diffs[4,2]+0.85,
              -seniors_amce_diffs[3,2]+0.85,
              0+0.85,
              -seniors_amce_diffs[2,2]+0.85,
              seniors_amce_diffs[1,2]+0.85,
              -smokers_amce_diffs[5,2]+0.85,
              -smokers_amce_diffs[4,2]+0.85,
              -smokers_amce_diffs[3,2]+0.85,
              0+0.85,
              -smokers_amce_diffs[2,2]+0.85,
              smokers_amce_diffs[1,2]+0.85,
              -women_amce_diffs[5,2]+0.85,
              -women_amce_diffs[4,2]+0.85,
              -women_amce_diffs[3,2]+0.85,
              0+0.85,
              -women_amce_diffs[2,2]+0.85,
              women_amce_diffs[1,2]+0.85
         ),
         x1=c(weight_amces[10,3],
              weight_amces[9,3],
              weight_amces[8,3],
              weight_amces[7,3],
              weight_amces[6,3],
              weight_amces[5,3],
              0,
              0,
              weight_amces[4,3],
              weight_amces[3,3],
              weight_amces[2,3],
              weight_amces[1,3],
              health_amces[10,3],
              health_amces[9,3],
              health_amces[8,3],
              health_amces[7,3],
              health_amces[6,3],
              health_amces[5,3],
              0,
              0,
              health_amces[4,3],
              health_amces[3,3],
              health_amces[2,3],
              health_amces[1,3],
              students_amces[10,3],
              students_amces[9,3],
              students_amces[8,3],
              students_amces[7,3],
              students_amces[6,3],
              students_amces[5,3],
              0,
              0,
              students_amces[4,3],
              students_amces[3,3],
              students_amces[2,3],
              students_amces[1,3],
              seniors_amces[10,3],
              seniors_amces[9,3],
              seniors_amces[8,3],
              seniors_amces[7,3],
              seniors_amces[6,3],
              seniors_amces[5,3],
              0,
              0,
              seniors_amces[4,3],
              seniors_amces[3,3],
              seniors_amces[2,3],
              seniors_amces[1,3],
              smokers_amces[10,3],
              smokers_amces[9,3],
              smokers_amces[8,3],
              smokers_amces[7,3],
              smokers_amces[6,3],
              smokers_amces[5,3],
              0,
              0,
              smokers_amces[4,3],
              smokers_amces[3,3],
              smokers_amces[2,3],
              smokers_amces[1,3],
              women_amces[10,3],
              women_amces[9,3],
              women_amces[8,3],
              women_amces[7,3],
              women_amces[6,3],
              women_amces[5,3],
              0,
              0,
              women_amces[4,3],
              women_amces[3,3],
              women_amces[2,3],
              women_amces[1,3],
              -weight_amce_diffs[5,3]+0.85,
              -weight_amce_diffs[4,3]+0.85,
              -weight_amce_diffs[3,3]+0.85,
              0+0.85,
              -weight_amce_diffs[2,3]+0.85,
              -weight_amce_diffs[1,3]+0.85,
              -health_amce_diffs[5,3]+0.85,
              -health_amce_diffs[4,3]+0.85,
              -health_amce_diffs[3,3]+0.85,
              0+0.85,
              -health_amce_diffs[2,3]+0.85,
              health_amce_diffs[1,3]+0.85,
              -students_amce_diffs[5,3]+0.85,
              -students_amce_diffs[4,3]+0.85,
              -students_amce_diffs[3,3]+0.85,
              0+0.85,
              -students_amce_diffs[2,3]+0.85,
              students_amce_diffs[1,3]+0.85,
              -seniors_amce_diffs[5,3]+0.85,
              seniors_amce_diffs[4,3]+0.85,
              -seniors_amce_diffs[3,3]+0.85,
              0+0.85,
              -seniors_amce_diffs[2,3]+0.85,
              seniors_amce_diffs[1,3]+0.85,
              -smokers_amce_diffs[5,3]+0.85,
              -smokers_amce_diffs[4,3]+0.85,
              -smokers_amce_diffs[3,3]+0.85,
              0+0.85,
              -smokers_amce_diffs[2,3]+0.85,
              smokers_amce_diffs[1,3]+0.85,
              -women_amce_diffs[5,3]+0.85,
              -women_amce_diffs[4,3]+0.85,
              -women_amce_diffs[3,3]+0.85,
              0+0.85,
              -women_amce_diffs[2,3]+0.85,
              women_amce_diffs[1,3]+0.85),
         y0=c(1.3, 1.7,
              2.3, 2.7,
              3.3, 3.7,
              4.3, 4.7,
              5.3, 5.7,
              6.3, 6.7,
              8.3, 8.7,
              9.3, 9.7,
              10.3, 10.7,
              11.3, 11.7,
              12.3, 12.7,
              13.3, 13.7,
              15.3, 15.7,
              16.3, 16.7,
              17.3, 17.7,
              18.3, 18.7,
              19.3, 19.7,
              20.3, 20.7,
              22.3, 22.7,
              23.3, 23.7,
              24.3, 24.7,
              25.3, 25.7,
              26.3, 26.7,
              27.3, 27.7,
              29.3, 29.7,
              30.3, 30.7,
              31.3, 31.7,
              32.3, 32.7,
              33.3, 33.7,
              34.3, 34.7,
              36.3, 36.7,
              37.3, 37.7,
              38.3, 38.7,
              39.3, 39.7,
              40.3, 40.7,
              41.3, 41.7,
              seq(1.5,6.5,1),
              seq(8.5,13.5,1),
              seq(15.5,20.5,1),
              seq(22.5,27.5,1),
              seq(29.5,34.5,1),
              seq(36.5,41.5,1)))
axis(1,at=c(-.25, 0, .25, .5), 
     labels = c("-25%", "0%", "25%", "50%"), cex.axis = 1.5)
axis(1,at=c(.60, .85, 1.10), 
     labels = c("-25%", "0%", "25%"), cex.axis = 1.5)
par(mgp=c(0,8,0))
axis(2,at=c(seq(1.5,6.5,1),
            seq(8.5,13.5,1),
            seq(15.5,20.5,1),
            seq(22.5,27.5,1),
            seq(29.5,34.5,1),
            seq(36.5,41.5,1)),
     las=2,
     labels=rep(c("Friendly-Relevant", "Friendly-Irrelevant",
                  "Neutral-Relevant", "Neutral-Irrelevant",
                  "Unfriendly-Relevant", "Unfriendly-Irrelevant"),6),
     tck=0,
     lwd = 0,
     line = 0,
     cex.axis=1, hadj=0)
brackets(x1=-.53, x2=-.53, y1=1.1, y2=6.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Losing Weight")), side = 2, at=4, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=8.1, y2=13.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Health Ins.")), side = 2, at=11, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=15.1, y2=20.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Students")), side = 2, at=18, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=22.1, y2=27.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Seniors")), side = 2, at=25, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=29.1, y2=34.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Smokers")), side = 2, at=32, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=36.1, y2=41.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Women")), side = 2, at=39, line = 9.2, cex = 1.75)
par(mgp=c(0,11,0))
mtext("Diff. from Baseline", side = 1, at=0.125, line = 4, cex = 1.75)
mtext("Magnitude of Under/Overestimation", side = 1, at=0.85, line = 4, cex = 1.75)
mtext("AMCEs by\n Outcome Options", side = 3, at=0.125, line=0.5, cex = 2)
mtext("Under/Overestimation with\nForced Choice", side = 3, at=0.85, line=0.5, cex = 2)
par(mar=c(0,0,0,0))
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
legend(x=0, y=0.8, legend=c("Forced Choice",
                            "Abstention Option",
                            "Under/Overestimation"), lty=1, lwd=1, 
       pch = c(19,17,15), col=c("black", "gray50", "black"), cex = 1.25)
dev.off()

############
# TABLE SM.8
############

# In order to make a table that resembles Table 4 in Mummolo (2016), we estimate
# separate models for the forced choice and abstention arms.

women_forced_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Forced"),])
summary(women_forced_reg)
women_forced_reg_vcovCL <- cluster.vcov(women_forced_reg, full_data_part$ResponseId[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Forced")])
coeftest(women_forced_reg, women_forced_reg_vcovCL)
women_forced_reg_ses <- coeftest(women_forced_reg, women_forced_reg_vcovCL)[,2]
women_forced_reg_pvals <- coeftest(women_forced_reg, women_forced_reg_vcovCL)[,4]
women_abstention_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Abstention"),])
summary(women_abstention_reg)
women_abstention_reg_vcovCL <- cluster.vcov(women_abstention_reg, full_data_part$ResponseId[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Abstention")])
coeftest(women_abstention_reg, women_abstention_reg_vcovCL)
women_abstention_reg_ses <- coeftest(women_abstention_reg, women_abstention_reg_vcovCL)[,2]
women_abstention_reg_pvals <- coeftest(women_abstention_reg, women_abstention_reg_vcovCL)[,4]

smokers_forced_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Forced"),])
summary(smokers_forced_reg)
smokers_forced_reg_vcovCL <- cluster.vcov(smokers_forced_reg, full_data_part$ResponseId[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Forced")])
coeftest(smokers_forced_reg, smokers_forced_reg_vcovCL)
smokers_forced_reg_ses <- coeftest(smokers_forced_reg, smokers_forced_reg_vcovCL)[,2]
smokers_forced_reg_pvals <- coeftest(smokers_forced_reg, smokers_forced_reg_vcovCL)[,4]
smokers_abstention_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Abstention"),])
summary(smokers_abstention_reg)
smokers_abstention_reg_vcovCL <- cluster.vcov(smokers_abstention_reg, full_data_part$ResponseId[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Abstention")])
coeftest(smokers_abstention_reg, smokers_abstention_reg_vcovCL)
smokers_abstention_reg_ses <- coeftest(smokers_abstention_reg, smokers_abstention_reg_vcovCL)[,2]
smokers_abstention_reg_pvals <- coeftest(smokers_abstention_reg, smokers_abstention_reg_vcovCL)[,4]

seniors_forced_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Forced"),])
summary(seniors_forced_reg)
seniors_forced_reg_vcovCL <- cluster.vcov(seniors_forced_reg, full_data_part$ResponseId[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Forced")])
coeftest(seniors_forced_reg, seniors_forced_reg_vcovCL)
seniors_forced_reg_ses <- coeftest(seniors_forced_reg, seniors_forced_reg_vcovCL)[,2]
seniors_forced_reg_pvals <- coeftest(seniors_forced_reg, seniors_forced_reg_vcovCL)[,4]
seniors_abstention_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Abstention"),])
summary(seniors_abstention_reg)
seniors_abstention_reg_vcovCL <- cluster.vcov(seniors_abstention_reg, full_data_part$ResponseId[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Abstention")])
coeftest(seniors_abstention_reg, seniors_abstention_reg_vcovCL)
seniors_abstention_reg_ses <- coeftest(seniors_abstention_reg, seniors_abstention_reg_vcovCL)[,2]
seniors_abstention_reg_pvals <- coeftest(seniors_abstention_reg, seniors_abstention_reg_vcovCL)[,4]

students_forced_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Forced"),])
summary(students_forced_reg)
students_forced_reg_vcovCL <- cluster.vcov(students_forced_reg, full_data_part$ResponseId[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Forced")])
coeftest(students_forced_reg, students_forced_reg_vcovCL)
students_forced_reg_ses <- coeftest(students_forced_reg, students_forced_reg_vcovCL)[,2]
students_forced_reg_pvals <- coeftest(students_forced_reg, students_forced_reg_vcovCL)[,4]
students_abstention_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Abstention"),])
summary(students_abstention_reg)
students_abstention_reg_vcovCL <- cluster.vcov(students_abstention_reg, full_data_part$ResponseId[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Abstention")])
coeftest(students_abstention_reg, students_abstention_reg_vcovCL)
students_abstention_reg_ses <- coeftest(students_abstention_reg, students_abstention_reg_vcovCL)[,2]
students_abstention_reg_pvals <- coeftest(students_abstention_reg, students_abstention_reg_vcovCL)[,4]

health_forced_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Forced"),])
summary(health_forced_reg)
health_forced_reg_vcovCL <- cluster.vcov(health_forced_reg, full_data_part$ResponseId[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Forced")])
coeftest(health_forced_reg, health_forced_reg_vcovCL)
health_forced_reg_ses <- coeftest(health_forced_reg, health_forced_reg_vcovCL)[,2]
health_forced_reg_pvals <- coeftest(health_forced_reg, health_forced_reg_vcovCL)[,4]
health_abstention_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Abstention"),])
summary(health_abstention_reg)
health_abstention_reg_vcovCL <- cluster.vcov(health_abstention_reg, full_data_part$ResponseId[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Abstention")])
coeftest(health_abstention_reg, health_abstention_reg_vcovCL)
health_abstention_reg_ses <- coeftest(health_abstention_reg, health_abstention_reg_vcovCL)[,2]
health_abstention_reg_pvals <- coeftest(health_abstention_reg, health_abstention_reg_vcovCL)[,4]

weight_forced_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Forced"),])
summary(weight_forced_reg)
weight_forced_reg_vcovCL <- cluster.vcov(weight_forced_reg, full_data_part$ResponseId[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Forced")])
coeftest(weight_forced_reg, weight_forced_reg_vcovCL)
weight_forced_reg_ses <- coeftest(weight_forced_reg, weight_forced_reg_vcovCL)[,2]
weight_forced_reg_pvals <- coeftest(weight_forced_reg, weight_forced_reg_vcovCL)[,4]
weight_abstention_reg <- lm(choice_bin ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Abstention"),])
summary(weight_abstention_reg)
weight_abstention_reg_vcovCL <- cluster.vcov(weight_abstention_reg, full_data_part$ResponseId[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Abstention")])
coeftest(weight_abstention_reg, weight_abstention_reg_vcovCL)
weight_abstention_reg_ses <- coeftest(weight_abstention_reg, weight_abstention_reg_vcovCL)[,2]
weight_abstention_reg_pvals <- coeftest(weight_abstention_reg, weight_abstention_reg_vcovCL)[,4]

print(
  texreg(l=list(women_forced_reg, women_abstention_reg,
              smokers_forced_reg, smokers_abstention_reg,
              seniors_forced_reg, seniors_abstention_reg,
              students_forced_reg, students_abstention_reg,
              health_forced_reg, health_abstention_reg,
              weight_forced_reg, weight_abstention_reg),
       custom.coef.names = c("Intercept", "Relevant Topic",
                             "Friendly Source", "Unfriendly Source",
                             "Relevant x Friendly", "Relevant x Unfriendly"),
       custom.model.names = c("Forced", "Abstention","Forced", "Abstention",
                              "Forced", "Abstention","Forced", "Abstention",
                              "Forced", "Abstention","Forced", "Abstention"),
       override.se = list(women_forced_reg_ses, women_abstention_reg_ses,
                          smokers_forced_reg_ses, smokers_abstention_reg_ses,
                          seniors_forced_reg_ses, seniors_abstention_reg_ses,
                          students_forced_reg_ses, students_abstention_reg_ses,
                          health_forced_reg_ses, health_abstention_reg_ses,
                          weight_forced_reg_ses, weight_abstention_reg_ses),
       override.pvalues = list(women_forced_reg_pvals, women_abstention_reg_pvals,
                               smokers_forced_reg_pvals, smokers_abstention_reg_pvals,
                               seniors_forced_reg_pvals, seniors_abstention_reg_pvals,
                               students_forced_reg_pvals, students_abstention_reg_pvals,
                               health_forced_reg_pvals, health_abstention_reg_pvals,
                               weight_forced_reg_pvals, weight_abstention_reg_pvals),
       caption = "Mummolo (2016) Replication Results",
       caption.above = TRUE,
       stars = 0.05,
       include.rsquared = FALSE,
       include.adjrs = FALSE)
  , file = "../tables/table_SM8.tex"
)
############
# TABLE SM.9
############

# Creating a new outcome variable for which abstentions are coded as NAs for
# both profiles in the task instead of zeroes
full_data_part$choice_bin_NA <- ifelse(full_data_part$choice=="Neither", NA, full_data_part$choice_bin)

# Repeating the same analysis as we conducted to make Table SM.8

women_forced_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Forced"),])
summary(women_forced_reg_NA)
women_forced_reg_NA_vcovCL <- cluster.vcov(women_forced_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Forced")])
coeftest(women_forced_reg_NA, women_forced_reg_NA_vcovCL)
women_forced_reg_NA_ses <- coeftest(women_forced_reg_NA, women_forced_reg_NA_vcovCL)[,2]
women_forced_reg_NA_pvals <- coeftest(women_forced_reg_NA, women_forced_reg_NA_vcovCL)[,4]
women_abstention_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Abstention"),])
summary(women_abstention_reg_NA)
women_abstention_reg_NA_vcovCL <- cluster.vcov(women_abstention_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_women==1 & full_data_part$exp_arm=="Abstention")])
coeftest(women_abstention_reg_NA, women_abstention_reg_NA_vcovCL)
women_abstention_reg_NA_ses <- coeftest(women_abstention_reg_NA, women_abstention_reg_NA_vcovCL)[,2]
women_abstention_reg_NA_pvals <- coeftest(women_abstention_reg_NA, women_abstention_reg_NA_vcovCL)[,4]

smokers_forced_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Forced"),])
summary(smokers_forced_reg_NA)
smokers_forced_reg_NA_vcovCL <- cluster.vcov(smokers_forced_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Forced")])
coeftest(smokers_forced_reg_NA, smokers_forced_reg_NA_vcovCL)
smokers_forced_reg_NA_ses <- coeftest(smokers_forced_reg_NA, smokers_forced_reg_NA_vcovCL)[,2]
smokers_forced_reg_NA_pvals <- coeftest(smokers_forced_reg_NA, smokers_forced_reg_NA_vcovCL)[,4]
smokers_abstention_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Abstention"),])
summary(smokers_abstention_reg_NA)
smokers_abstention_reg_NA_vcovCL <- cluster.vcov(smokers_abstention_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_smokers==1 & full_data_part$exp_arm=="Abstention")])
coeftest(smokers_abstention_reg_NA, smokers_abstention_reg_NA_vcovCL)
smokers_abstention_reg_NA_ses <- coeftest(smokers_abstention_reg_NA, smokers_abstention_reg_NA_vcovCL)[,2]
smokers_abstention_reg_NA_pvals <- coeftest(smokers_abstention_reg_NA, smokers_abstention_reg_NA_vcovCL)[,4]

seniors_forced_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Forced"),])
summary(seniors_forced_reg_NA)
seniors_forced_reg_NA_vcovCL <- cluster.vcov(seniors_forced_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Forced")])
coeftest(seniors_forced_reg_NA, seniors_forced_reg_NA_vcovCL)
seniors_forced_reg_NA_ses <- coeftest(seniors_forced_reg_NA, seniors_forced_reg_NA_vcovCL)[,2]
seniors_forced_reg_NA_pvals <- coeftest(seniors_forced_reg_NA, seniors_forced_reg_NA_vcovCL)[,4]
seniors_abstention_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Abstention"),])
summary(seniors_abstention_reg_NA)
seniors_abstention_reg_NA_vcovCL <- cluster.vcov(seniors_abstention_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_seniors==1 & full_data_part$exp_arm=="Abstention")])
coeftest(seniors_abstention_reg_NA, seniors_abstention_reg_NA_vcovCL)
seniors_abstention_reg_NA_ses <- coeftest(seniors_abstention_reg_NA, seniors_abstention_reg_NA_vcovCL)[,2]
seniors_abstention_reg_NA_pvals <- coeftest(seniors_abstention_reg_NA, seniors_abstention_reg_NA_vcovCL)[,4]

students_forced_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Forced"),])
summary(students_forced_reg_NA)
students_forced_reg_NA_vcovCL <- cluster.vcov(students_forced_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Forced")])
coeftest(students_forced_reg_NA, students_forced_reg_NA_vcovCL)
students_forced_reg_NA_ses <- coeftest(students_forced_reg_NA, students_forced_reg_NA_vcovCL)[,2]
students_forced_reg_NA_pvals <- coeftest(students_forced_reg_NA, students_forced_reg_NA_vcovCL)[,4]
students_abstention_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Abstention"),])
summary(students_abstention_reg_NA)
students_abstention_reg_NA_vcovCL <- cluster.vcov(students_abstention_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_students==1 & full_data_part$exp_arm=="Abstention")])
coeftest(students_abstention_reg_NA, students_abstention_reg_NA_vcovCL)
students_abstention_reg_NA_ses <- coeftest(students_abstention_reg_NA, students_abstention_reg_NA_vcovCL)[,2]
students_abstention_reg_NA_pvals <- coeftest(students_abstention_reg_NA, students_abstention_reg_NA_vcovCL)[,4]

health_forced_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Forced"),])
summary(health_forced_reg_NA)
health_forced_reg_NA_vcovCL <- cluster.vcov(health_forced_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Forced")])
coeftest(health_forced_reg_NA, health_forced_reg_NA_vcovCL)
health_forced_reg_NA_ses <- coeftest(health_forced_reg_NA, health_forced_reg_NA_vcovCL)[,2]
health_forced_reg_NA_pvals <- coeftest(health_forced_reg_NA, health_forced_reg_NA_vcovCL)[,4]
health_abstention_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Abstention"),])
summary(health_abstention_reg_NA)
health_abstention_reg_NA_vcovCL <- cluster.vcov(health_abstention_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_health==1 & full_data_part$exp_arm=="Abstention")])
coeftest(health_abstention_reg_NA, health_abstention_reg_NA_vcovCL)
health_abstention_reg_NA_ses <- coeftest(health_abstention_reg_NA, health_abstention_reg_NA_vcovCL)[,2]
health_abstention_reg_NA_pvals <- coeftest(health_abstention_reg_NA, health_abstention_reg_NA_vcovCL)[,4]

weight_forced_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Forced"),])
summary(weight_forced_reg_NA)
weight_forced_reg_NA_vcovCL <- cluster.vcov(weight_forced_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Forced")])
coeftest(weight_forced_reg_NA, weight_forced_reg_NA_vcovCL)
weight_forced_reg_NA_ses <- coeftest(weight_forced_reg_NA, weight_forced_reg_NA_vcovCL)[,2]
weight_forced_reg_NA_pvals <- coeftest(weight_forced_reg_NA, weight_forced_reg_NA_vcovCL)[,4]
weight_abstention_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness, data = full_data_part[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Abstention"),])
summary(weight_abstention_reg_NA)
weight_abstention_reg_NA_vcovCL <- cluster.vcov(weight_abstention_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_weight==1 & full_data_part$exp_arm=="Abstention")])
coeftest(weight_abstention_reg_NA, weight_abstention_reg_NA_vcovCL)
weight_abstention_reg_NA_ses <- coeftest(weight_abstention_reg_NA, weight_abstention_reg_NA_vcovCL)[,2]
weight_abstention_reg_NA_pvals <- coeftest(weight_abstention_reg_NA, weight_abstention_reg_NA_vcovCL)[,4]

print(
  texreg(l=list(women_forced_reg_NA, women_abstention_reg_NA,
              smokers_forced_reg_NA, smokers_abstention_reg_NA,
              seniors_forced_reg_NA, seniors_abstention_reg_NA,
              students_forced_reg_NA, students_abstention_reg_NA,
              health_forced_reg_NA, health_abstention_reg_NA,
              weight_forced_reg_NA, weight_abstention_reg_NA),
       custom.coef.names = c("Intercept", "Relevant Topic",
                             "Friendly Source", "Unfriendly Source",
                             "Relevant x Friendly", "Relevant x Unfriendly"),
       custom.model.names = c("Forced", "Abstention","Forced", "Abstention",
                              "Forced", "Abstention","Forced", "Abstention",
                              "Forced", "Abstention","Forced", "Abstention"),
       override.se = list(women_forced_reg_NA_ses, women_abstention_reg_NA_ses,
                          smokers_forced_reg_NA_ses, smokers_abstention_reg_NA_ses,
                          seniors_forced_reg_NA_ses, seniors_abstention_reg_NA_ses,
                          students_forced_reg_NA_ses, students_abstention_reg_NA_ses,
                          health_forced_reg_NA_ses, health_abstention_reg_NA_ses,
                          weight_forced_reg_NA_ses, weight_abstention_reg_NA_ses),
       override.pvalues = list(women_forced_reg_NA_pvals, women_abstention_reg_NA_pvals,
                               smokers_forced_reg_NA_pvals, smokers_abstention_reg_NA_pvals,
                               seniors_forced_reg_NA_pvals, seniors_abstention_reg_NA_pvals,
                               students_forced_reg_NA_pvals, students_abstention_reg_NA_pvals,
                               health_forced_reg_NA_pvals, health_abstention_reg_NA_pvals,
                               weight_forced_reg_NA_pvals, weight_abstention_reg_NA_pvals),
       caption = "Mummolo (2016) Replication Results",
       caption.above = TRUE,
       stars = 0.05,
       include.rsquared = FALSE,
       include.adjrs = FALSE)
  , file = "../tables/table_SM9.tex"
)

#############
# FIGURE SM.6
#############

# The below code is identical to that used to create Figure SM.5 above EXCEPT
# it uses the "abstention as NA" outcome measure

women_int_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_women==1),])
summary(women_int_reg_NA)
women_int_reg_NA_vcovCL <- cluster.vcov(women_int_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_women==1)])
coeftest(women_int_reg_NA, women_int_reg_NA_vcovCL)
women_int_reg_NA <-withVCov(women_int_reg_NA, vcov = women_int_reg_NA_vcovCL)
women_int_reg_NA_ses <- coeftest(women_int_reg_NA, women_int_reg_NA_vcovCL)[,2]
women_int_reg_NA_pvals <- coeftest(women_int_reg_NA, women_int_reg_NA_vcovCL)[,4]

smokers_int_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_smokers==1),])
summary(smokers_int_reg_NA)
smokers_int_reg_NA_vcovCL <- cluster.vcov(smokers_int_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_smokers==1)])
coeftest(smokers_int_reg_NA, smokers_int_reg_NA_vcovCL)
smokers_int_reg_NA <-withVCov(smokers_int_reg_NA, vcov = smokers_int_reg_NA_vcovCL)
smokers_int_reg_NA_ses <- coeftest(smokers_int_reg_NA, smokers_int_reg_NA_vcovCL)[,2]
smokers_int_reg_NA_pvals <- coeftest(smokers_int_reg_NA, smokers_int_reg_NA_vcovCL)[,4]

seniors_int_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_seniors==1),])
summary(seniors_int_reg_NA)
seniors_int_reg_NA_vcovCL <- cluster.vcov(seniors_int_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_seniors==1)])
coeftest(seniors_int_reg_NA, seniors_int_reg_NA_vcovCL)
seniors_int_reg_NA <-withVCov(seniors_int_reg_NA, vcov = seniors_int_reg_NA_vcovCL)
seniors_int_reg_NA_ses <- coeftest(seniors_int_reg_NA, seniors_int_reg_NA_vcovCL)[,2]
seniors_int_reg_NA_pvals <- coeftest(seniors_int_reg_NA, seniors_int_reg_NA_vcovCL)[,4]

students_int_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_students==1),])
summary(students_int_reg_NA)
students_int_reg_NA_vcovCL <- cluster.vcov(students_int_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_students==1)])
coeftest(students_int_reg_NA, students_int_reg_NA_vcovCL)
students_int_reg_NA <-withVCov(students_int_reg_NA, vcov = students_int_reg_NA_vcovCL)
students_int_reg_NA_ses <- coeftest(students_int_reg_NA, students_int_reg_NA_vcovCL)[,2]
students_int_reg_NA_pvals <- coeftest(students_int_reg_NA, students_int_reg_NA_vcovCL)[,4]

health_int_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_health==1),])
summary(health_int_reg_NA)
health_int_reg_NA_vcovCL <- cluster.vcov(health_int_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_health==1)])
coeftest(health_int_reg_NA, health_int_reg_NA_vcovCL)
health_int_reg_NA <-withVCov(health_int_reg_NA, vcov = health_int_reg_NA_vcovCL)
health_int_reg_NA_ses <- coeftest(health_int_reg_NA, health_int_reg_NA_vcovCL)[,2]
health_int_reg_NA_pvals <- coeftest(health_int_reg_NA, health_int_reg_NA_vcovCL)[,4]

weight_int_reg_NA <- lm(choice_bin_NA ~ topic_relevance*source_friendliness*exp_arm, data = full_data_part[which(full_data_part$ap_weight==1),])
summary(weight_int_reg_NA)
weight_int_reg_NA_vcovCL <- cluster.vcov(weight_int_reg_NA, full_data_part$ResponseId[which(full_data_part$ap_weight==1)])
coeftest(weight_int_reg_NA, weight_int_reg_NA_vcovCL)
weight_int_reg_NA <-withVCov(weight_int_reg_NA, vcov = weight_int_reg_NA_vcovCL)
weight_int_reg_NA_ses <- coeftest(weight_int_reg_NA, weight_int_reg_NA_vcovCL)[,2]
weight_int_reg_NA_pvals <- coeftest(weight_int_reg_NA, weight_int_reg_NA_vcovCL)[,4]

women_amces_NA <- confint(glht(women_int_reg_NA, 
                               linfct = c("source_friendlinessUnfriendly = 0",
                                          "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                          "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                          "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                          "topic_relevanceRelevant = 0",
                                          "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                          "source_friendlinessFriendly = 0",
                                          "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                          "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                          "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                          calpha = univariate_calpha())$confint

women_amce_NA_diffs <- confint(glht(women_int_reg_NA, 
                                    linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                               "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                               "topic_relevanceRelevant:exp_armAbstention = 0",
                                               "source_friendlinessFriendly:exp_armAbstention = 0",
                                               "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                               calpha = 2.94)$confint

smokers_amces_NA <- confint(glht(smokers_int_reg_NA, 
                                 linfct = c("source_friendlinessUnfriendly = 0",
                                            "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                            "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                            "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                            "topic_relevanceRelevant = 0",
                                            "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                            "source_friendlinessFriendly = 0",
                                            "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                            "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                            "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                            calpha = univariate_calpha())$confint

smokers_amce_NA_diffs <- confint(glht(smokers_int_reg_NA, 
                                      linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                 "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                 "topic_relevanceRelevant:exp_armAbstention = 0",
                                                 "source_friendlinessFriendly:exp_armAbstention = 0",
                                                 "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                                 calpha = 2.94)$confint

seniors_amces_NA <- confint(glht(seniors_int_reg_NA, 
                                 linfct = c("source_friendlinessUnfriendly = 0",
                                            "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                            "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                            "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                            "topic_relevanceRelevant = 0",
                                            "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                            "source_friendlinessFriendly = 0",
                                            "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                            "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                            "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                            calpha = univariate_calpha())$confint

seniors_amce_NA_diffs <- confint(glht(seniors_int_reg_NA, 
                                      linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                 "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                 "topic_relevanceRelevant:exp_armAbstention = 0",
                                                 "source_friendlinessFriendly:exp_armAbstention = 0",
                                                 "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                                 calpha = 2.94)$confint

students_amces_NA <- confint(glht(students_int_reg_NA, 
                                  linfct = c("source_friendlinessUnfriendly = 0",
                                             "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                             "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                             "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                             "topic_relevanceRelevant = 0",
                                             "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                             "source_friendlinessFriendly = 0",
                                             "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                             "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                             "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                             calpha = univariate_calpha())$confint

students_amce_NA_diffs <- confint(glht(students_int_reg_NA, 
                                       linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                  "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                  "topic_relevanceRelevant:exp_armAbstention = 0",
                                                  "source_friendlinessFriendly:exp_armAbstention = 0",
                                                  "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                                  calpha = 2.94)$confint

health_amces_NA <- confint(glht(health_int_reg_NA, 
                                linfct = c("source_friendlinessUnfriendly = 0",
                                           "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                           "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                           "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                           "topic_relevanceRelevant = 0",
                                           "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                           "source_friendlinessFriendly = 0",
                                           "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                           "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                           "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                           calpha = univariate_calpha())$confint

health_amce_NA_diffs <- confint(glht(health_int_reg_NA, 
                                     linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                "topic_relevanceRelevant:exp_armAbstention = 0",
                                                "source_friendlinessFriendly:exp_armAbstention = 0",
                                                "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                                calpha = 2.94)$confint

weight_amces_NA <- confint(glht(weight_int_reg_NA, 
                                linfct = c("source_friendlinessUnfriendly = 0",
                                           "source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention = 0",
                                           "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly = 0",
                                           "source_friendlinessUnfriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessUnfriendly + source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                           "topic_relevanceRelevant = 0",
                                           "topic_relevanceRelevant + topic_relevanceRelevant:exp_armAbstention = 0",
                                           "source_friendlinessFriendly = 0",
                                           "source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention = 0",
                                           "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly = 0",
                                           "source_friendlinessFriendly + topic_relevanceRelevant + topic_relevanceRelevant:source_friendlinessFriendly + source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                           calpha = univariate_calpha())$confint

weight_amce_NA_diffs <- confint(glht(weight_int_reg_NA, 
                                     linfct = c("source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                "source_friendlinessUnfriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessUnfriendly:exp_armAbstention = 0",
                                                "topic_relevanceRelevant:exp_armAbstention = 0",
                                                "source_friendlinessFriendly:exp_armAbstention = 0",
                                                "source_friendlinessFriendly:exp_armAbstention + topic_relevanceRelevant:exp_armAbstention + topic_relevanceRelevant:source_friendlinessFriendly:exp_armAbstention = 0")),
                                calpha = 2.94)$confint

pdf(file = "../figures/figure_SM6.pdf", family = "Times", height = 18, width=13)
layout(matrix(c(1,2), nrow=2, byrow=TRUE), heights = c(0.85,0.15))
par(mar=c(6.1,12,5.75,.5))
plot(x=c(weight_amces_NA[10,1],
         weight_amces_NA[9,1],
         weight_amces_NA[8,1],
         weight_amces_NA[7,1],
         weight_amces_NA[6,1],
         weight_amces_NA[5,1],
         0,
         0,
         weight_amces_NA[4,1],
         weight_amces_NA[3,1],
         weight_amces_NA[2,1],
         weight_amces_NA[1,1],
         health_amces_NA[10,1],
         health_amces_NA[9,1],
         health_amces_NA[8,1],
         health_amces_NA[7,1],
         health_amces_NA[6,1],
         health_amces_NA[5,1],
         0,
         0,
         health_amces_NA[4,1],
         health_amces_NA[3,1],
         health_amces_NA[2,1],
         health_amces_NA[1,1],
         students_amces_NA[10,1],
         students_amces_NA[9,1],
         students_amces_NA[8,1],
         students_amces_NA[7,1],
         students_amces_NA[6,1],
         students_amces_NA[5,1],
         0,
         0,
         students_amces_NA[4,1],
         students_amces_NA[3,1],
         students_amces_NA[2,1],
         students_amces_NA[1,1],
         seniors_amces_NA[10,1],
         seniors_amces_NA[9,1],
         seniors_amces_NA[8,1],
         seniors_amces_NA[7,1],
         seniors_amces_NA[6,1],
         seniors_amces_NA[5,1],
         0,
         0,
         seniors_amces_NA[4,1],
         seniors_amces_NA[3,1],
         seniors_amces_NA[2,1],
         seniors_amces_NA[1,1],
         smokers_amces_NA[10,1],
         smokers_amces_NA[9,1],
         smokers_amces_NA[8,1],
         smokers_amces_NA[7,1],
         smokers_amces_NA[6,1],
         smokers_amces_NA[5,1],
         0,
         0,
         smokers_amces_NA[4,1],
         smokers_amces_NA[3,1],
         smokers_amces_NA[2,1],
         smokers_amces_NA[1,1],
         women_amces_NA[10,1],
         women_amces_NA[9,1],
         women_amces_NA[8,1],
         women_amces_NA[7,1],
         women_amces_NA[6,1],
         women_amces_NA[5,1],
         0,
         0,
         women_amces_NA[4,1],
         women_amces_NA[3,1],
         women_amces_NA[2,1],
         women_amces_NA[1,1],
         -weight_amce_NA_diffs[5,1]+0.85,
         -weight_amce_NA_diffs[4,1]+0.85,
         -weight_amce_NA_diffs[3,1]+0.85,
         0+0.85,
         -weight_amce_NA_diffs[2,1]+0.85,
         -weight_amce_NA_diffs[1,1]+0.85,
         -health_amce_NA_diffs[5,1]+0.85,
         -health_amce_NA_diffs[4,1]+0.85,
         -health_amce_NA_diffs[3,1]+0.85,
         0+0.85,
         -health_amce_NA_diffs[2,1]+0.85,
         health_amce_NA_diffs[1,1]+0.85,
         -students_amce_NA_diffs[5,1]+0.85,
         -students_amce_NA_diffs[4,1]+0.85,
         -students_amce_NA_diffs[3,1]+0.85,
         0+0.85,
         -students_amce_NA_diffs[2,1]+0.85,
         students_amce_NA_diffs[1,1]+0.85,
         -seniors_amce_NA_diffs[5,1]+0.85,
         seniors_amce_NA_diffs[4,1]+0.85,
         -seniors_amce_NA_diffs[3,1]+0.85,
         0+0.85,
         -seniors_amce_NA_diffs[2,1]+0.85,
         seniors_amce_NA_diffs[1,1]+0.85,
         -smokers_amce_NA_diffs[5,1]+0.85,
         -smokers_amce_NA_diffs[4,1]+0.85,
         -smokers_amce_NA_diffs[3,1]+0.85,
         0+0.85,
         -smokers_amce_NA_diffs[2,1]+0.85,
         smokers_amce_NA_diffs[1,1]+0.85,
         -women_amce_NA_diffs[5,1]+0.85,
         -women_amce_NA_diffs[4,1]+0.85,
         -women_amce_NA_diffs[3,1]+0.85,
         0+0.85,
         -women_amce_NA_diffs[2,1]+0.85,
         women_amce_NA_diffs[1,1]+0.85), 
     y=c(1.3, 1.7,
         2.3, 2.7,
         3.3, 3.7,
         4.3, 4.7,
         5.3, 5.7,
         6.3, 6.7,
         8.3, 8.7,
         9.3, 9.7,
         10.3, 10.7,
         11.3, 11.7,
         12.3, 12.7,
         13.3, 13.7,
         15.3, 15.7,
         16.3, 16.7,
         17.3, 17.7,
         18.3, 18.7,
         19.3, 19.7,
         20.3, 20.7,
         22.3, 22.7,
         23.3, 23.7,
         24.3, 24.7,
         25.3, 25.7,
         26.3, 26.7,
         27.3, 27.7,
         29.3, 29.7,
         30.3, 30.7,
         31.3, 31.7,
         32.3, 32.7,
         33.3, 33.7,
         34.3, 34.7,
         36.3, 36.7,
         37.3, 37.7,
         38.3, 38.7,
         39.3, 39.7,
         40.3, 40.7,
         41.3, 41.7,
         seq(1.5,6.5,1),
         seq(8.5,13.5,1),
         seq(15.5,20.5,1),
         seq(22.5,27.5,1),
         seq(29.5,34.5,1),
         seq(36.5,41.5,1)),
     xlim=c(-0.25,1.15),
     xaxt="n",
     ylim=c(1, 42),
     pch=c(rep(c(17,19),36), rep(15, 36)),
     col=c(rep(c("gray50","black"),36), rep("black", 36)),
     tck=-.02,
     cex.axis=0.9,
     cex=1.25,
     ylab="",
     yaxt="n",
     xlab="",
     axes = FALSE,
     panel.first = c(abline(v=0,lwd=2, col="gray70",lty=2),
                     abline(v=0.85,lwd=2, col="gray70",lty=2)))
segments(col=c(rep(c("gray50", "black"),36), rep("black", 36)),
         x0=c(weight_amces_NA[10,2],
              weight_amces_NA[9,2],
              weight_amces_NA[8,2],
              weight_amces_NA[7,2],
              weight_amces_NA[6,2],
              weight_amces_NA[5,2],
              0,
              0,
              weight_amces_NA[4,2],
              weight_amces_NA[3,2],
              weight_amces_NA[2,2],
              weight_amces_NA[1,2],
              health_amces_NA[10,2],
              health_amces_NA[9,2],
              health_amces_NA[8,2],
              health_amces_NA[7,2],
              health_amces_NA[6,2],
              health_amces_NA[5,2],
              0,
              0,
              health_amces_NA[4,2],
              health_amces_NA[3,2],
              health_amces_NA[2,2],
              health_amces_NA[1,2],
              students_amces_NA[10,2],
              students_amces_NA[9,2],
              students_amces_NA[8,2],
              students_amces_NA[7,2],
              students_amces_NA[6,2],
              students_amces_NA[5,2],
              0,
              0,
              students_amces_NA[4,2],
              students_amces_NA[3,2],
              students_amces_NA[2,2],
              students_amces_NA[1,2],
              seniors_amces_NA[10,2],
              seniors_amces_NA[9,2],
              seniors_amces_NA[8,2],
              seniors_amces_NA[7,2],
              seniors_amces_NA[6,2],
              seniors_amces_NA[5,2],
              0,
              0,
              seniors_amces_NA[4,2],
              seniors_amces_NA[3,2],
              seniors_amces_NA[2,2],
              seniors_amces_NA[1,2],
              smokers_amces_NA[10,2],
              smokers_amces_NA[9,2],
              smokers_amces_NA[8,2],
              smokers_amces_NA[7,2],
              smokers_amces_NA[6,2],
              smokers_amces_NA[5,2],
              0,
              0,
              smokers_amces_NA[4,2],
              smokers_amces_NA[3,2],
              smokers_amces_NA[2,2],
              smokers_amces_NA[1,2],
              women_amces_NA[10,2],
              women_amces_NA[9,2],
              women_amces_NA[8,2],
              women_amces_NA[7,2],
              women_amces_NA[6,2],
              women_amces_NA[5,2],
              0,
              0,
              women_amces_NA[4,2],
              women_amces_NA[3,2],
              women_amces_NA[2,2],
              women_amces_NA[1,2],
              -weight_amce_NA_diffs[5,2]+0.85,
              -weight_amce_NA_diffs[4,2]+0.85,
              -weight_amce_NA_diffs[3,2]+0.85,
              0+0.85,
              -weight_amce_NA_diffs[2,2]+0.85,
              -weight_amce_NA_diffs[1,2]+0.85,
              -health_amce_NA_diffs[5,2]+0.85,
              -health_amce_NA_diffs[4,2]+0.85,
              -health_amce_NA_diffs[3,2]+0.85,
              0+0.85,
              -health_amce_NA_diffs[2,2]+0.85,
              health_amce_NA_diffs[1,2]+0.85,
              -students_amce_NA_diffs[5,2]+0.85,
              -students_amce_NA_diffs[4,2]+0.85,
              -students_amce_NA_diffs[3,2]+0.85,
              0+0.85,
              -students_amce_NA_diffs[2,2]+0.85,
              students_amce_NA_diffs[1,2]+0.85,
              -seniors_amce_NA_diffs[5,2]+0.85,
              seniors_amce_NA_diffs[4,2]+0.85,
              -seniors_amce_NA_diffs[3,2]+0.85,
              0+0.85,
              -seniors_amce_NA_diffs[2,2]+0.85,
              seniors_amce_NA_diffs[1,2]+0.85,
              -smokers_amce_NA_diffs[5,2]+0.85,
              -smokers_amce_NA_diffs[4,2]+0.85,
              -smokers_amce_NA_diffs[3,2]+0.85,
              0+0.85,
              -smokers_amce_NA_diffs[2,2]+0.85,
              smokers_amce_NA_diffs[1,2]+0.85,
              -women_amce_NA_diffs[5,2]+0.85,
              -women_amce_NA_diffs[4,2]+0.85,
              -women_amce_NA_diffs[3,2]+0.85,
              0+0.85,
              -women_amce_NA_diffs[2,2]+0.85,
              women_amce_NA_diffs[1,2]+0.85
         ),
         x1=c(weight_amces_NA[10,3],
              weight_amces_NA[9,3],
              weight_amces_NA[8,3],
              weight_amces_NA[7,3],
              weight_amces_NA[6,3],
              weight_amces_NA[5,3],
              0,
              0,
              weight_amces_NA[4,3],
              weight_amces_NA[3,3],
              weight_amces_NA[2,3],
              weight_amces_NA[1,3],
              health_amces_NA[10,3],
              health_amces_NA[9,3],
              health_amces_NA[8,3],
              health_amces_NA[7,3],
              health_amces_NA[6,3],
              health_amces_NA[5,3],
              0,
              0,
              health_amces_NA[4,3],
              health_amces_NA[3,3],
              health_amces_NA[2,3],
              health_amces_NA[1,3],
              students_amces_NA[10,3],
              students_amces_NA[9,3],
              students_amces_NA[8,3],
              students_amces_NA[7,3],
              students_amces_NA[6,3],
              students_amces_NA[5,3],
              0,
              0,
              students_amces_NA[4,3],
              students_amces_NA[3,3],
              students_amces_NA[2,3],
              students_amces_NA[1,3],
              seniors_amces_NA[10,3],
              seniors_amces_NA[9,3],
              seniors_amces_NA[8,3],
              seniors_amces_NA[7,3],
              seniors_amces_NA[6,3],
              seniors_amces_NA[5,3],
              0,
              0,
              seniors_amces_NA[4,3],
              seniors_amces_NA[3,3],
              seniors_amces_NA[2,3],
              seniors_amces_NA[1,3],
              smokers_amces_NA[10,3],
              smokers_amces_NA[9,3],
              smokers_amces_NA[8,3],
              smokers_amces_NA[7,3],
              smokers_amces_NA[6,3],
              smokers_amces_NA[5,3],
              0,
              0,
              smokers_amces_NA[4,3],
              smokers_amces_NA[3,3],
              smokers_amces_NA[2,3],
              smokers_amces_NA[1,3],
              women_amces_NA[10,3],
              women_amces_NA[9,3],
              women_amces_NA[8,3],
              women_amces_NA[7,3],
              women_amces_NA[6,3],
              women_amces_NA[5,3],
              0,
              0,
              women_amces_NA[4,3],
              women_amces_NA[3,3],
              women_amces_NA[2,3],
              women_amces_NA[1,3],
              -weight_amce_NA_diffs[5,3]+0.85,
              -weight_amce_NA_diffs[4,3]+0.85,
              -weight_amce_NA_diffs[3,3]+0.85,
              0+0.85,
              -weight_amce_NA_diffs[2,3]+0.85,
              -weight_amce_NA_diffs[1,3]+0.85,
              -health_amce_NA_diffs[5,3]+0.85,
              -health_amce_NA_diffs[4,3]+0.85,
              -health_amce_NA_diffs[3,3]+0.85,
              0+0.85,
              -health_amce_NA_diffs[2,3]+0.85,
              health_amce_NA_diffs[1,3]+0.85,
              -students_amce_NA_diffs[5,3]+0.85,
              -students_amce_NA_diffs[4,3]+0.85,
              -students_amce_NA_diffs[3,3]+0.85,
              0+0.85,
              -students_amce_NA_diffs[2,3]+0.85,
              students_amce_NA_diffs[1,3]+0.85,
              -seniors_amce_NA_diffs[5,3]+0.85,
              seniors_amce_NA_diffs[4,3]+0.85,
              -seniors_amce_NA_diffs[3,3]+0.85,
              0+0.85,
              -seniors_amce_NA_diffs[2,3]+0.85,
              seniors_amce_NA_diffs[1,3]+0.85,
              -smokers_amce_NA_diffs[5,3]+0.85,
              -smokers_amce_NA_diffs[4,3]+0.85,
              -smokers_amce_NA_diffs[3,3]+0.85,
              0+0.85,
              -smokers_amce_NA_diffs[2,3]+0.85,
              smokers_amce_NA_diffs[1,3]+0.85,
              -women_amce_NA_diffs[5,3]+0.85,
              -women_amce_NA_diffs[4,3]+0.85,
              -women_amce_NA_diffs[3,3]+0.85,
              0+0.85,
              -women_amce_NA_diffs[2,3]+0.85,
              women_amce_NA_diffs[1,3]+0.85),
         y0=c(1.3, 1.7,
              2.3, 2.7,
              3.3, 3.7,
              4.3, 4.7,
              5.3, 5.7,
              6.3, 6.7,
              8.3, 8.7,
              9.3, 9.7,
              10.3, 10.7,
              11.3, 11.7,
              12.3, 12.7,
              13.3, 13.7,
              15.3, 15.7,
              16.3, 16.7,
              17.3, 17.7,
              18.3, 18.7,
              19.3, 19.7,
              20.3, 20.7,
              22.3, 22.7,
              23.3, 23.7,
              24.3, 24.7,
              25.3, 25.7,
              26.3, 26.7,
              27.3, 27.7,
              29.3, 29.7,
              30.3, 30.7,
              31.3, 31.7,
              32.3, 32.7,
              33.3, 33.7,
              34.3, 34.7,
              36.3, 36.7,
              37.3, 37.7,
              38.3, 38.7,
              39.3, 39.7,
              40.3, 40.7,
              41.3, 41.7,
              seq(1.5,6.5,1),
              seq(8.5,13.5,1),
              seq(15.5,20.5,1),
              seq(22.5,27.5,1),
              seq(29.5,34.5,1),
              seq(36.5,41.5,1)))
axis(1,at=c(-.25, 0, .25, .5), 
     labels = c("-25%", "0%", "25%", "50%"), cex.axis = 1.5)
axis(1,at=c(.60, .85, 1.10), 
     labels = c("-25%", "0%", "25%"), cex.axis = 1.5)
par(mgp=c(0,8,0))
axis(2,at=c(seq(1.5,6.5,1),
            seq(8.5,13.5,1),
            seq(15.5,20.5,1),
            seq(22.5,27.5,1),
            seq(29.5,34.5,1),
            seq(36.5,41.5,1)),
     las=2,
     labels=rep(c("Friendly-Relevant", "Friendly-Irrelevant",
                  "Neutral-Relevant", "Neutral-Irrelevant",
                  "Unfriendly-Relevant", "Unfriendly-Irrelevant"),6),
     tck=0,
     lwd = 0,
     line = 0,
     cex.axis=1, hadj=0)
brackets(x1=-.53, x2=-.53, y1=1.1, y2=6.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Losing Weight")), side = 2, at=4, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=8.1, y2=13.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Health Ins.")), side = 2, at=11, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=15.1, y2=20.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Students")), side = 2, at=18, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=22.1, y2=27.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Seniors")), side = 2, at=25, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=29.1, y2=34.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Smokers")), side = 2, at=32, line = 9.2, cex = 1.75)
brackets(x1=-.53, x2=-.53, y1=36.1, y2=41.9, type=4, ticks = NA, xpd=TRUE, h=0.02)
mtext(expression(~bold("Women")), side = 2, at=39, line = 9.2, cex = 1.75)
par(mgp=c(0,11,0))
mtext("Diff. from Baseline", side = 1, at=0.125, line = 4, cex = 1.75)
mtext("Magnitude of Under/Overestimation", side = 1, at=0.85, line = 4, cex = 1.75)
mtext("AMCEs by\n Outcome Options", side = 3, at=0.125, line=0.5, cex = 2)
mtext("Under/Overestimation with\nForced Choice", side = 3, at=0.85, line=0.5, cex = 2)
par(mar=c(0,0,0,0))
plot(0,0, type="n", axes=FALSE, xlab="", ylab="")
legend(x=0, y=0.8, legend=c("Forced Choice",
                            "Abstention Option",
                            "Under/Overestimation"), lty=1, lwd=1, 
       pch = c(19,17,15), col=c("black", "gray50", "black"), cex = 1.25)
dev.off()

#############
# FIGURE SM.7
#############

# The following code chunks take the information about both profiles in each task
# and code task-level information about whether one or both stories in the
# task were relevant to the respondent and from a friendly source, then
# binds the task-level dataframes into a single long dataframe

long_data1 <- subset(data, select = c(ResponseId, headline1_A, source1_A,
                                      headline1_B, source1_B, choice1))
long_data1$task <- 1
colnames(long_data1) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")

long_data2 <- subset(data, select = c(ResponseId, headline2_A, source2_A,
                                      headline2_B, source2_B, choice2))
long_data2$task <- 2
colnames(long_data2) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")

long_data3 <- subset(data, select = c(ResponseId, headline3_A, source3_A,
                                      headline3_B, source3_B, choice3))
long_data3$task <- 3
colnames(long_data3) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")

long_data4 <- subset(data, select = c(ResponseId, headline4_A, source4_A,
                                      headline4_B, source4_B, choice4))
long_data4$task <- 4
colnames(long_data4) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")

long_data5 <- subset(data, select = c(ResponseId, headline5_A, source5_A,
                                      headline5_B, source5_B, choice5))
long_data5$task <- 5
colnames(long_data5) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")

long_data6 <- subset(data, select = c(ResponseId, headline6_A, source6_A,
                                      headline6_B, source6_B, choice6))
long_data6$task <- 6
colnames(long_data6) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")

long_data7 <- subset(data, select = c(ResponseId, headline7_A, source7_A,
                                      headline7_B, source7_B, choice7))
long_data7$task <- 7
colnames(long_data7) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")


long_data8 <- subset(data, select = c(ResponseId, headline8_A, source8_A,
                                      headline8_B, source8_B, choice8))
long_data8$task <- 8
colnames(long_data8) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")


long_data9 <- subset(data, select = c(ResponseId, headline9_A, source9_A,
                                      headline9_B, source9_B, choice9))
long_data9$task <- 9
colnames(long_data9) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                          "sourceB", "choice", "task")


long_data10 <- subset(data, select = c(ResponseId, headline10_A, source10_A,
                                       headline10_B, source10_B, choice10))
long_data10$task <- 10
colnames(long_data10) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                           "sourceB", "choice", "task")


long_data11 <- subset(data, select = c(ResponseId, headline11_A, source11_A,
                                       headline11_B, source11_B, choice11))
long_data11$task <- 11
colnames(long_data11) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                           "sourceB", "choice", "task")


long_data12 <- subset(data, select = c(ResponseId, headline12_A, source12_A,
                                       headline12_B, source12_B, choice12))
long_data12$task <- 12
colnames(long_data12) <- c("ResponseId", "headlineA", "sourceA", "headlineB",
                           "sourceB", "choice", "task")

long_data <- rbind(long_data1, long_data2, long_data3, long_data4,
                   long_data5, long_data6, long_data7, long_data8,
                   long_data9, long_data10, long_data11, long_data12)

# Merge back in respondent demographic characteristics
full_data_task <- merge(demo_data, long_data, by = "ResponseId")

# The following code chunks ascertain for each profile the friendliness of the
# source and the relevant of each topic
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="MSNBC" & full_data_task$pid7>4, "Friendly", NA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="MSNBC" & full_data_task$pid7<4, "Unfriendly", full_data_task$source_friendlinessA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="MSNBC" & full_data_task$pid7==4, "Neutral", full_data_task$source_friendlinessA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="Fox News" & full_data_task$pid7<4, "Friendly", full_data_task$source_friendlinessA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="Fox News" & full_data_task$pid7>4, "Unfriendly", full_data_task$source_friendlinessA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="Fox News" & full_data_task$pid7==4, "Neutral", full_data_task$source_friendlinessA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="USA Today" & full_data_task$pid7<4, "Neutral", full_data_task$source_friendlinessA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="USA Today" & full_data_task$pid7>4, "Neutral", full_data_task$source_friendlinessA)
full_data_task$source_friendlinessA <- ifelse(full_data_task$sourceA=="USA Today" & full_data_task$pid7==4, "Neutral", full_data_task$source_friendlinessA)

full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="2021 MLB Baseball Season Guide", "Irrelevant", NA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="2021 MLB Baseball Season Outlook", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="2021 MLB Baseball Season Preview", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Abortion restrictions reinstated by federal court" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Abortion restrictions reinstated by federal court" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Bill to guarantee equal pay for women blocked in Senate" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Bill to guarantee equal pay for women blocked in Senate" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Bill to mandate equal pay for women fails in Senate vote" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Bill to mandate equal pay for women fails in Senate vote" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Celebrities' worst dating stories", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Celebrity dating fails", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Celebrity dating nightmares", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Federal court puts abortion restrictions back in place" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Federal court puts abortion restrictions back in place" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Federal court reinstates abortion restrictions" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Federal court reinstates abortion restrictions" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Congress considering cuts to Social Security" & full_data_task$ap_seniors==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Congress considering cuts to Social Security" & full_data_task$ap_seniors==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Congress debates cuts to Social Security" & full_data_task$ap_seniors==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Congress debates cuts to Social Security" & full_data_task$ap_seniors==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Congress weighs cuts to Social Security" & full_data_task$ap_seniors==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Congress weighs cuts to Social Security" & full_data_task$ap_seniors==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="How to drop all that student loan debt" & full_data_task$ap_students==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="How to drop all that student loan debt" & full_data_task$ap_students==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Here are the new members of Joe Biden's Cabinet", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Introducing Joe Biden's new Cabinet members", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Meet the new members of Joe Biden's Cabinet", "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Obamacare enrollment figures met with doubt" & full_data_task$ap_health==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Obamacare enrollment figures met with doubt" & full_data_task$ap_health==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Obamacare enrollment numbers called into question" & full_data_task$ap_health==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Obamacare enrollment numbers called into question" & full_data_task$ap_health==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Some question Obamacare enrollment figures" & full_data_task$ap_health==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Some question Obamacare enrollment figures" & full_data_task$ap_health==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Senate votes against bill that would ensure equal pay for women" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Senate votes against bill that would ensure equal pay for women" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Smokers who quit may cut heart risk faster than had been thought, study finds" & full_data_task$ap_smokers==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Smokers who quit may cut heart risk faster than had been thought, study finds" & full_data_task$ap_smokers==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study finds potential cause of breast cancer" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study finds potential cause of breast cancer" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study finds those who quit smoking repair heart damage fast than expected" & full_data_task$ap_smokers==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study finds those who quit smoking repair heart damage fast than expected" & full_data_task$ap_smokers==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study identifies possible breast cancer cause" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study identifies possible breast cancer cause" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study points to new cause of breast cancer" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study points to new cause of breast cancer" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study: quitting smoking repairs heart damage more quickly than was thought" & full_data_task$ap_smokers==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Study: quitting smoking repairs heart damage more quickly than was thought" & full_data_task$ap_smokers==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="The best ways to get rid of student loan debt" & full_data_task$ap_students==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="The best ways to get rid of student loan debt" & full_data_task$ap_students==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Tips for getting rid of student loan debt" & full_data_task$ap_students==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Tips for getting rid of student loan debt" & full_data_task$ap_students==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Weight-loss tips and success stories" & full_data_task$ap_weight==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Weight-loss tips and success stories" & full_data_task$ap_weight==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Weight-loss tips that make a difference" & full_data_task$ap_weight==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Weight-loss tips that make a difference" & full_data_task$ap_weight==0, "Irrelevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Weight-loss tips that work" & full_data_task$ap_weight==1, "Relevant", full_data_task$topic_relevanceA)
full_data_task$topic_relevanceA <- ifelse(full_data_task$headlineA=="Weight-loss tips that work" & full_data_task$ap_weight==0, "Irrelevant", full_data_task$topic_relevanceA)

full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="MSNBC" & full_data_task$pid7>4, "Friendly", NA)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="MSNBC" & full_data_task$pid7<4, "Unfriendly", full_data_task$source_friendlinessB)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="MSNBC" & full_data_task$pid7==4, "Neutral", full_data_task$source_friendlinessB)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="Fox News" & full_data_task$pid7<4, "Friendly", full_data_task$source_friendlinessB)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="Fox News" & full_data_task$pid7>4, "Unfriendly", full_data_task$source_friendlinessB)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="Fox News" & full_data_task$pid7==4, "Neutral", full_data_task$source_friendlinessB)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="USA Today" & full_data_task$pid7<4, "Neutral", full_data_task$source_friendlinessB)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="USA Today" & full_data_task$pid7>4, "Neutral", full_data_task$source_friendlinessB)
full_data_task$source_friendlinessB <- ifelse(full_data_task$sourceB=="USA Today" & full_data_task$pid7==4, "Neutral", full_data_task$source_friendlinessB)

full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="2021 MLB Baseball Season Guide", "Irrelevant", NA)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="2021 MLB Baseball Season Outlook", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="2021 MLB Baseball Season Preview", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Abortion restrictions reinstated by federal court" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Abortion restrictions reinstated by federal court" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Bill to guarantee equal pay for women blocked in Senate" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Bill to guarantee equal pay for women blocked in Senate" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Bill to mandate equal pay for women fails in Senate vote" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Bill to mandate equal pay for women fails in Senate vote" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Celebrities' worst dating stories", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Celebrity dating fails", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Celebrity dating nightmares", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Federal court puts abortion restrictions back in place" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Federal court puts abortion restrictions back in place" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Federal court reinstates abortion restrictions" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Federal court reinstates abortion restrictions" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Congress considering cuts to Social Security" & full_data_task$ap_seniors==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Congress considering cuts to Social Security" & full_data_task$ap_seniors==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Congress debates cuts to Social Security" & full_data_task$ap_seniors==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Congress debates cuts to Social Security" & full_data_task$ap_seniors==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Congress weighs cuts to Social Security" & full_data_task$ap_seniors==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Congress weighs cuts to Social Security" & full_data_task$ap_seniors==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="How to drop all that student loan debt" & full_data_task$ap_students==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="How to drop all that student loan debt" & full_data_task$ap_students==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Here are the new members of Joe Biden's Cabinet", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Introducing Joe Biden's new Cabinet members", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Meet the new members of Joe Biden's Cabinet", "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Obamacare enrollment figures met with doubt" & full_data_task$ap_health==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Obamacare enrollment figures met with doubt" & full_data_task$ap_health==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Obamacare enrollment numbers called into question" & full_data_task$ap_health==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Obamacare enrollment numbers called into question" & full_data_task$ap_health==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Some question Obamacare enrollment figures" & full_data_task$ap_health==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Some question Obamacare enrollment figures" & full_data_task$ap_health==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Senate votes against bill that would ensure equal pay for women" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Senate votes against bill that would ensure equal pay for women" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Smokers who quit may cut heart risk faster than had been thought, study finds" & full_data_task$ap_smokers==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Smokers who quit may cut heart risk faster than had been thought, study finds" & full_data_task$ap_smokers==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study finds potential cause of breast cancer" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study finds potential cause of breast cancer" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study finds those who quit smoking repair heart damage fast than expected" & full_data_task$ap_smokers==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study finds those who quit smoking repair heart damage fast than expected" & full_data_task$ap_smokers==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study identifies possible breast cancer cause" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study identifies possible breast cancer cause" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study points to new cause of breast cancer" & full_data_task$ap_women==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study points to new cause of breast cancer" & full_data_task$ap_women==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study: quitting smoking repairs heart damage more quickly than was thought" & full_data_task$ap_smokers==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Study: quitting smoking repairs heart damage more quickly than was thought" & full_data_task$ap_smokers==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="The best ways to get rid of student loan debt" & full_data_task$ap_students==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="The best ways to get rid of student loan debt" & full_data_task$ap_students==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Tips for getting rid of student loan debt" & full_data_task$ap_students==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Tips for getting rid of student loan debt" & full_data_task$ap_students==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Weight-loss tips and success stories" & full_data_task$ap_weight==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Weight-loss tips and success stories" & full_data_task$ap_weight==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Weight-loss tips that make a difference" & full_data_task$ap_weight==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Weight-loss tips that make a difference" & full_data_task$ap_weight==0, "Irrelevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Weight-loss tips that work" & full_data_task$ap_weight==1, "Relevant", full_data_task$topic_relevanceB)
full_data_task$topic_relevanceB <- ifelse(full_data_task$headlineB=="Weight-loss tips that work" & full_data_task$ap_weight==0, "Irrelevant", full_data_task$topic_relevanceB)

# Coding whether in the profile each respondent chose a story or abstained
full_data_task$abstention <- ifelse(full_data_task$choice=="News Selection A" | full_data_task$choice=="News Selection B", 0, NA)
full_data_task$abstention <- ifelse(full_data_task$choice=="Neither", 1, full_data_task$abstention)

# Coding for each task whether both stories were friendly, both unfriendly,
# or whether one was friendly and the other was unfriendly
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Friendly" & full_data_task$source_friendlinessB=="Friendly",
                                "Both Friendly", NA)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Neutral" & full_data_task$source_friendlinessB=="Neutral",
                                "Both Neutral", full_data_task$source_pair)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Unfriendly" & full_data_task$source_friendlinessB=="Unfriendly",
                                "Both Unfriendly", full_data_task$source_pair)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Friendly" & full_data_task$source_friendlinessB=="Neutral",
                                "Friendly-Neutral", full_data_task$source_pair)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Neutral" & full_data_task$source_friendlinessB=="Friendly",
                                "Friendly-Neutral", full_data_task$source_pair)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Friendly" & full_data_task$source_friendlinessB=="Unfriendly",
                                "Friendly-Unfriendly", full_data_task$source_pair)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Unfriendly" & full_data_task$source_friendlinessB=="Friendly",
                                "Friendly-Unfriendly", full_data_task$source_pair)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Neutral" & full_data_task$source_friendlinessB=="Unfriendly",
                                "Neutral-Unfriendly", full_data_task$source_pair)
full_data_task$source_pair <- ifelse(full_data_task$source_friendlinessA=="Unfriendly" & full_data_task$source_friendlinessB=="Neutral",
                                "Neutral-Unfriendly", full_data_task$source_pair)

# Coding for each task whether both stories were relevant, whether both were
# irrelevant, or whether one was relevant and the other was irrelevant
full_data_task$topic_pair <- ifelse(full_data_task$topic_relevanceA=="Relevant" & full_data_task$topic_relevanceB=="Relevant",
                               "Both Relevant", NA)
full_data_task$topic_pair <- ifelse(full_data_task$topic_relevanceA=="Relevant" & full_data_task$topic_relevanceB=="Irrelevant",
                               "Relevant-Irrelevant", full_data_task$topic_pair)
full_data_task$topic_pair <- ifelse(full_data_task$topic_relevanceA=="Irrelevant" & full_data_task$topic_relevanceB=="Relevant",
                               "Relevant-Irrelevant", full_data_task$topic_pair)
full_data_task$topic_pair <- ifelse(full_data_task$topic_relevanceA=="Irrelevant" & full_data_task$topic_relevanceB=="Irrelevant",
                               "Both Irrelevant", full_data_task$topic_pair)

# Mummolo retained only partisans, so all who identified as true independents
# need to be dropped
full_data_task_part <- full_data_task[which(full_data_task$pid7!=4),]

# Specifying levels and order of levels for factor variables
full_data_task_part$source_pair <- factor(full_data_task_part$source_pair, levels = c("Both Unfriendly",
                                                                            "Neutral-Unfriendly",
                                                                            "Both Neutral",
                                                                            "Friendly-Unfriendly",
                                                                            "Friendly-Neutral",
                                                                            "Both Friendly"))
full_data_task_part$topic_pair <- factor(full_data_task_part$topic_pair, levels = c("Both Irrelevant",
                                                                          "Relevant-Irrelevant",
                                                                          "Both Relevant"))

# Regressing abstention on the task-level characteritics
why_abstention <- lm(abstention ~ source_pair*topic_pair, data = full_data_task_part[which(full_data_task_part$exp_arm=="Abstention"),])
summary(why_abstention)

# Calculating the predicted probability of abstention for the specified
# task-level characteristics
why_abstention_predictions <- predict(why_abstention, newdata = data.frame("source_pair"=rep(c("Both Unfriendly",
                                                                                               "Neutral-Unfriendly",
                                                                                               "Both Neutral",
                                                                                               "Friendly-Unfriendly",
                                                                                               "Friendly-Neutral",
                                                                                               "Both Friendly"
                                                                              ),3),
                                                                              "topic_pair"=c(rep("Both Irrelevant",6),
                                                                                             rep("Relevant-Irrelevant",6),
                                                                                             rep("Both Relevant",6))),
                                                                              se.fit = TRUE)

pdf(file = "../figures/figure_SM7.pdf", family = "Times", height = 8, width=8)
par(mar=c(6.1,14,2,1))
plot(x=c(why_abstention_predictions$fit[1:18]), 
     y=c(1:6, 8:13, 15:20),
     xlim=c(0,0.50),
     xaxt="n",
     ylim=c(1, 20),
     pch=19,
     col="black",
     tck=-.02,
     cex.axis=0.9,
     cex=1.25,
     ylab="",
     yaxt="n",
     xlab="",
     axes = FALSE)
segments(col="black",
         x0=c(why_abstention_predictions$fit[1:18]-1.96*why_abstention_predictions$se.fit[1:18]),
         x1=c(why_abstention_predictions$fit[1:18]+1.96*why_abstention_predictions$se.fit[1:18]),
         y0=c(1:6, 8:13, 15:20))
axis(1,at=c(0.00, 0.10, 0.20, 0.30, 0.40, 0.50), 
     labels = c("0.00", "0.10", "0.20", "0.30", "0.40", "0.50"), cex.axis = 1.5)
par(mgp=c(0,8,0))
axis(2,at=c(1:6, 8:13, 15:20),
     las=2,
     labels=rep(c("Both Unfriendly",
                  "Neutral-Unfriendly",
                  "Both Neutral",
                  "Friendly-Unfriendly",
                  "Friendly-Neutral",
                  "Both Friendly"),3),
     tck=0,
     lwd = 0,
     line = 1,
     cex.axis=1, hadj=0)
brackets(x1=-.02, x2=-.02, y1=0.51, y2=6.49, type=4, ticks = NA, xpd=TRUE, h=0.01)
mtext(expression(~bold("Both Irrelevant")), side = 2, at=3.5, line = 10, cex = 1.25)
brackets(x1=-.02, x2=-.02, y1=7.51, y2=13.49, type=4, ticks = NA, xpd=TRUE, h=0.01)
mtext(expression(~bold("Relevant-Irrelevant")), side = 2, at=10.5, line = 10, cex = 1.25)
brackets(x1=-.02, x2=-.02, y1=14.51, y2=20.49, type=4, ticks = NA, xpd=TRUE, h=0.01)
mtext(expression(~bold("Both Relevant")), side = 2, at=17.5, line = 10, cex = 1.25)
par(mgp=c(0,11,0))
mtext("Probability of Abstention", side = 1, at=0.25, line = 4, cex = 1.75)
dev.off()

#############
# TABLE SM.10
#############

# in the respondent-profile dataframe, creating a binary indicator for
# abstention
full_data_task_part$abstention <- ifelse(full_data_task_part$choice=="News Selection A" | full_data_task_part$choice=="News Selection B", 0, NA)
full_data_task_part$abstention <- ifelse(full_data_task_part$choice=="Neither", 1, full_data_task_part$abstention)

# cleaning demographic data for each respondent, subsetting to a demographics
# dataframe, then merging with the respondent-profile dataframe

data$age <- 2021 - as.numeric(data$Q2)

data$female <- ifelse(data$Q4=="Female", 1, NA)
data$female <- ifelse(data$Q4=="Male", 0, data$female)
data$female <- ifelse(data$Q4=="Other (Please specify)", 0, data$female)

data$education <- ifelse(data$Q6=="Some high school, or less", 1, NA)
data$education <- ifelse(data$Q6=="High school graduate or GED", 2, data$education)
data$education <- ifelse(data$Q6=="Some college, no 4-year degree", 3, data$education)
data$education <- ifelse(data$Q6=="College graduate", 4, data$education)
data$education <- ifelse(data$Q6=="Post-graduate degree", 5, data$education)

data$black <- ifelse(data$Q7=="Black or African-American", 1, NA)
data$black <- ifelse(data$Q7!="Black or African-American" & data$Q7!="", 0, data$black)

data$asian <- ifelse(data$Q7=="Asian", 1, NA)
data$asian <- ifelse(data$Q7!="Asian" & data$Q7!="", 0, data$asian)

data$otherrace <- ifelse(data$Q7=="American Indian or Alaska Native", 1, NA)
data$otherrace <- ifelse(data$Q7=="Native Hawaiian or other Pacific Islander", 1, data$otherrace)
data$otherrace <- ifelse(data$Q7=="Other (Please specify)", 1, data$otherrace)
data$otherrace <- ifelse(data$Q7!="Other (Please specify)" & 
                           data$Q7!="Native Hawaiian or other Pacific Islander" & 
                           data$Q7!="American Indian or Alaska Native" & 
                           data$Q7!="", 0, data$otherrace)

data$hispanic <- ifelse(data$Q8=="Yes", 1, NA)
data$hispanic <- ifelse(data$Q8=="No", 0, data$hispanic)

data$income <- ifelse(data$Q9=="Less than $25,000", 1, NA)
data$income <- ifelse(data$Q9=="$25,000-$50,000", 2, data$income)
data$income <- ifelse(data$Q9=="$50,000-$75,000", 3, data$income)
data$income <- ifelse(data$Q9=="$75,000-$100,000", 4, data$income)
data$income <- ifelse(data$Q9=="$100,000-$200,000", 5, data$income)
data$income <- ifelse(data$Q9=="$200,000 or more", 6, data$income)

data$ideology <- ifelse(data$Q11=="Very conservative", 1, NA)
data$ideology <- ifelse(data$Q11=="Conservative", 2, data$ideology)
data$ideology <- ifelse(data$Q11=="Slightly conservative", 3, data$ideology)
data$ideology <- ifelse(data$Q11=="Moderate", 4, data$ideology)
data$ideology <- ifelse(data$Q11=="Slightly liberal", 5, data$ideology)
data$ideology <- ifelse(data$Q11=="Liberal", 6, data$ideology)
data$ideology <- ifelse(data$Q11=="Very liberal", 7, data$ideology)

data$pid3 <- ifelse(data$Q14=="Strong Republican", "Republican", NA)
data$pid3 <- ifelse(data$Q14=="Not a very strong Republican", "Republican", data$pid3)
data$pid3 <- ifelse(data$Q15=="Republican", "Republican", data$pid3)
data$pid3 <- ifelse(data$Q15=="Neither", "Independent", data$pid3)
data$pid3 <- ifelse(data$Q15=="Democratic", "Democrat", data$pid3)
data$pid3 <- ifelse(data$Q13=="Not a very strong Democrat", "Democrat", data$pid3)
data$pid3 <- ifelse(data$Q13=="Strong Democrat", "Democrat", data$pid3)
data$pid3 <- ifelse(data$Q12=="Other", "Other", data$pid3)

demo_data <- subset(data, select = c(ResponseId, age, female, 
                                     education, black, asian, otherrace, 
                                     hispanic, income, ideology))

full_data_task_part <- merge(demo_data, full_data_task_part, by = "ResponseId")

why_abstention_demo <- lm(abstention ~ education + black + asian + otherrace +
                       hispanic + income + ideology, data = full_data_task_part[which(full_data_task_part$exp_arm=="Abstention"),])
summary(why_abstention_demo)

print(
  texreg(l=list(why_abstention_demo),
       custom.coef.map = list("(Intercept)" = "Intercept",
                              "education" = "Education",
                              "asian" = "Asian",
                              "black" = "Black",
                              "otherrace" = "Other Non-White Race",
                              "hispanic" = "Hispanic",
                              "income" = "Income",
                              "ideology" = "Ideology"
       ),
       caption = "Demographic correlates of abstention in Mummolo replication",
       caption.above = TRUE,
       stars = 0.05,
       include.rsquared = FALSE,
       include.adjrs = FALSE)
  , file = "../tables/table_SM10.tex"
)
#################
# IN-TEXT RESULTS
#################

# PAGE 10, FOOTNOTE 7, AND SM.40--HOW COMMON ARE ABSTENTIONS?

# across all tasks 
table(full_data_task_part$abstention[which(full_data_task_part$exp_arm=="Abstention")])
prop.table(table(full_data_task_part$abstention[which(full_data_task_part$exp_arm=="Abstention")]))

# how many respondents abstained at least once?
ever_abstain <- aggregate(full_data_task_part$abstention[!is.na(full_data_task_part$abstention)], 
                          by = list(full_data_task_part$ResponseId[!is.na(full_data_task_part$abstention)]), FUN = "sum")
length(ever_abstain$Group.1[which(ever_abstain$x!=0)])
length(ever_abstain$Group.1[which(ever_abstain$x==0)])
length(ever_abstain$Group.1[which(ever_abstain$x!=0)])/dim(ever_abstain)[1]

# SM.44--HOW OFTEN DO RESPONDENTS SWITCH FROM CHOOSING TO ABSTAINING WHEN
# ALLOWED TO DO SO?

# Depending on whether respondents were in the forced choice or abstention arms,
# their choices in each task were stored as answers to different questions
# in Qualtrics.  Because their final 13th task can have been any of the
# 12 tasks they saw previously, we determine if they "switched" on their last
# task by examining the correspondence between their experimental arm and which
# of the two possible question slots for the task they had responses recorded;
# if for any task they had answers in both slots, that is the task for which
# they had the opportunity to switch

data$choice_initial1 <- ifelse(data$exp_arm=="Forced", data$Q23, data$Q52)
data$choice_final1 <- ifelse(data$exp_arm=="Abstention", data$Q23, data$Q52)

data$choice_initial2 <- ifelse(data$exp_arm=="Forced", data$Q25, data$Q54)
data$choice_final2 <- ifelse(data$exp_arm=="Abstention", data$Q25, data$Q54)

data$choice_initial3 <- ifelse(data$exp_arm=="Forced", data$Q29, data$Q56)
data$choice_final3 <- ifelse(data$exp_arm=="Abstention", data$Q29, data$Q56)

data$choice_initial4 <- ifelse(data$exp_arm=="Forced", data$Q31, data$Q58)
data$choice_final4 <- ifelse(data$exp_arm=="Abstention", data$Q31, data$Q58)

data$choice_initial5 <- ifelse(data$exp_arm=="Forced", data$Q33, data$Q60)
data$choice_final5 <- ifelse(data$exp_arm=="Abstention", data$Q33, data$Q60)

data$choice_initial6 <- ifelse(data$exp_arm=="Forced", data$Q35, data$Q62)
data$choice_final6 <- ifelse(data$exp_arm=="Abstention", data$Q35, data$Q62)

data$choice_initial7 <- ifelse(data$exp_arm=="Forced", data$Q37, data$Q64)
data$choice_final7 <- ifelse(data$exp_arm=="Abstention", data$Q37, data$Q64)

data$choice_initial8 <- ifelse(data$exp_arm=="Forced", data$Q39, data$Q66)
data$choice_final8 <- ifelse(data$exp_arm=="Abstention", data$Q39, data$Q66)

data$choice_initial9 <- ifelse(data$exp_arm=="Forced", data$Q41, data$Q68)
data$choice_final9 <- ifelse(data$exp_arm=="Abstention", data$Q41, data$Q68)

data$choice_initial10 <- ifelse(data$exp_arm=="Forced", data$Q43, data$Q70)
data$choice_final10 <- ifelse(data$exp_arm=="Abstention", data$Q43, data$Q70)

data$choice_initial11 <- ifelse(data$exp_arm=="Forced", data$Q45, data$Q72)
data$choice_final11 <- ifelse(data$exp_arm=="Abstention", data$Q45, data$Q72)

data$choice_initial12 <- ifelse(data$exp_arm=="Forced", data$Q27, data$Q74)
data$choice_final12 <- ifelse(data$exp_arm=="Abstention", data$Q27, data$Q74)

data$repeated_task <- ifelse(data$choice_initial1!="" & data$choice_final1!="", 1, NA)
data$repeated_task <- ifelse(data$choice_initial2!="" & data$choice_final2!="", 2, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial3!="" & data$choice_final3!="", 3, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial4!="" & data$choice_final4!="", 4, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial5!="" & data$choice_final5!="", 5, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial6!="" & data$choice_final6!="", 6, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial7!="" & data$choice_final7!="", 7, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial8!="" & data$choice_final8!="", 8, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial9!="" & data$choice_final9!="", 9, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial10!="" & data$choice_final10!="", 10, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial11!="" & data$choice_final11!="", 11, data$repeated_task)
data$repeated_task <- ifelse(data$choice_initial11!="" & data$choice_final12!="", 12, data$repeated_task)

data$initial_choice <- ifelse(data$repeated_task==1, data$choice_initial1, NA)
data$initial_choice <- ifelse(data$repeated_task==2, data$choice_initial2, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==3, data$choice_initial3, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==4, data$choice_initial4, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==5, data$choice_initial5, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==6, data$choice_initial6, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==7, data$choice_initial7, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==8, data$choice_initial8, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==9, data$choice_initial9, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==10, data$choice_initial10, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==11, data$choice_initial11, data$initial_choice)
data$initial_choice <- ifelse(data$repeated_task==12, data$choice_initial12, data$initial_choice)

data$final_choice <- ifelse(data$repeated_task==1, data$choice_final1, NA)
data$final_choice <- ifelse(data$repeated_task==2, data$choice_final2, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==3, data$choice_final3, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==4, data$choice_final4, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==5, data$choice_final5, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==6, data$choice_final6, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==7, data$choice_final7, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==8, data$choice_final8, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==9, data$choice_final9, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==10, data$choice_final10, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==11, data$choice_final11, data$final_choice)
data$final_choice <- ifelse(data$repeated_task==12, data$choice_final12, data$final_choice)

table(data$initial_choice[which(data$exp_arm=="Forced")], data$final_choice[which(data$exp_arm=="Forced")])
prop.table(table(data$initial_choice[which(data$exp_arm=="Forced")], data$final_choice[which(data$exp_arm=="Forced")]))
